Actual source code: rvector.c
petsc-3.11.4 2019-09-28
2: /*
3: Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
4: These are the vector functions the user calls.
5: */
6: #include <petsc/private/vecimpl.h>
7: static PetscInt VecGetSubVectorSavedStateId = -1;
9: PETSC_EXTERN PetscErrorCode VecValidValues(Vec vec,PetscInt argnum,PetscBool begin)
10: {
11: #if defined(PETSC_USE_DEBUG)
12: PetscErrorCode ierr;
13: PetscInt n,i;
14: const PetscScalar *x;
17: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL)
18: if ((vec->petscnative || vec->ops->getarray) && (vec->valid_GPU_array == PETSC_OFFLOAD_CPU || vec->valid_GPU_array == PETSC_OFFLOAD_BOTH)) {
19: #else
20: if (vec->petscnative || vec->ops->getarray) {
21: #endif
22: VecGetLocalSize(vec,&n);
23: VecGetArrayRead(vec,&x);
24: for (i=0; i<n; i++) {
25: if (begin) {
26: if (PetscIsInfOrNanScalar(x[i])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FP,"Vec entry at local location %D is not-a-number or infinite at beginning of function: Parameter number %D",i,argnum);
27: } else {
28: if (PetscIsInfOrNanScalar(x[i])) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FP,"Vec entry at local location %D is not-a-number or infinite at end of function: Parameter number %D",i,argnum);
29: }
30: }
31: VecRestoreArrayRead(vec,&x);
32: }
33: return(0);
34: #else
35: return 0;
36: #endif
37: }
39: /*@
40: VecMaxPointwiseDivide - Computes the maximum of the componentwise division max = max_i abs(x_i/y_i).
42: Logically Collective on Vec
44: Input Parameters:
45: . x, y - the vectors
47: Output Parameter:
48: . max - the result
50: Level: advanced
52: Notes:
53: x and y may be the same vector
54: if a particular y_i is zero, it is treated as 1 in the above formula
56: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs()
57: @*/
58: PetscErrorCode VecMaxPointwiseDivide(Vec x,Vec y,PetscReal *max)
59: {
69: VecCheckSameSize(x,1,y,2);
70: (*x->ops->maxpointwisedivide)(x,y,max);
71: return(0);
72: }
74: /*@
75: VecDot - Computes the vector dot product.
77: Collective on Vec
79: Input Parameters:
80: . x, y - the vectors
82: Output Parameter:
83: . val - the dot product
85: Performance Issues:
86: $ per-processor memory bandwidth
87: $ interprocessor latency
88: $ work load inbalance that causes certain processes to arrive much earlier than others
90: Notes for Users of Complex Numbers:
91: For complex vectors, VecDot() computes
92: $ val = (x,y) = y^H x,
93: where y^H denotes the conjugate transpose of y. Note that this corresponds to the usual "mathematicians" complex
94: inner product where the SECOND argument gets the complex conjugate. Since the BLASdot() complex conjugates the first
95: first argument we call the BLASdot() with the arguments reversed.
97: Use VecTDot() for the indefinite form
98: $ val = (x,y) = y^T x,
99: where y^T denotes the transpose of y.
101: Level: intermediate
103: Concepts: inner product
104: Concepts: vector^inner product
106: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDotRealPart()
107: @*/
108: PetscErrorCode VecDot(Vec x,Vec y,PetscScalar *val)
109: {
119: VecCheckSameSize(x,1,y,2);
121: PetscLogEventBegin(VEC_Dot,x,y,0,0);
122: (*x->ops->dot)(x,y,val);
123: PetscLogEventEnd(VEC_Dot,x,y,0,0);
124: return(0);
125: }
127: /*@
128: VecDotRealPart - Computes the real part of the vector dot product.
130: Collective on Vec
132: Input Parameters:
133: . x, y - the vectors
135: Output Parameter:
136: . val - the real part of the dot product;
138: Performance Issues:
139: $ per-processor memory bandwidth
140: $ interprocessor latency
141: $ work load inbalance that causes certain processes to arrive much earlier than others
143: Notes for Users of Complex Numbers:
144: See VecDot() for more details on the definition of the dot product for complex numbers
146: For real numbers this returns the same value as VecDot()
148: For complex numbers in C^n (that is a vector of n components with a complex number for each component) this is equal to the usual real dot product on the
149: the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)
151: Developer Note: This is not currently optimized to compute only the real part of the dot product.
153: Level: intermediate
155: Concepts: inner product
156: Concepts: vector^inner product
158: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecDot(), VecDotNorm2()
159: @*/
160: PetscErrorCode VecDotRealPart(Vec x,Vec y,PetscReal *val)
161: {
163: PetscScalar fdot;
166: VecDot(x,y,&fdot);
167: *val = PetscRealPart(fdot);
168: return(0);
169: }
171: /*@
172: VecNorm - Computes the vector norm.
174: Collective on Vec
176: Input Parameters:
177: + x - the vector
178: - type - one of NORM_1, NORM_2, NORM_INFINITY. Also available
179: NORM_1_AND_2, which computes both norms and stores them
180: in a two element array.
182: Output Parameter:
183: . val - the norm
185: Notes:
186: $ NORM_1 denotes sum_i |x_i|
187: $ NORM_2 denotes sqrt(sum_i |x_i|^2)
188: $ NORM_INFINITY denotes max_i |x_i|
190: For complex numbers NORM_1 will return the traditional 1 norm of the 2 norm of the complex numbers; that is the 1
191: norm of the absolutely values of the complex entries. In PETSc 3.6 and earlier releases it returned the 1 norm of
192: the 1 norm of the complex entries (what is returned by the BLAS routine asum()). Both are valid norms but most
193: people expect the former.
195: Level: intermediate
197: Performance Issues:
198: $ per-processor memory bandwidth
199: $ interprocessor latency
200: $ work load inbalance that causes certain processes to arrive much earlier than others
202: Concepts: norm
203: Concepts: vector^norm
205: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNormAvailable(),
206: VecNormBegin(), VecNormEnd()
208: @*/
209: PetscErrorCode VecNorm(Vec x,NormType type,PetscReal *val)
210: {
211: PetscBool flg;
219: /*
220: * Cached data?
221: */
222: if (type!=NORM_1_AND_2) {
223: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,flg);
224: if (flg) return(0);
225: }
226: PetscLogEventBegin(VEC_Norm,x,0,0,0);
227: (*x->ops->norm)(x,type,val);
228: PetscLogEventEnd(VEC_Norm,x,0,0,0);
230: if (type!=NORM_1_AND_2) {
231: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[type],*val);
232: }
233: return(0);
234: }
236: /*@
237: VecNormAvailable - Returns the vector norm if it is already known.
239: Not Collective
241: Input Parameters:
242: + x - the vector
243: - type - one of NORM_1, NORM_2, NORM_INFINITY. Also available
244: NORM_1_AND_2, which computes both norms and stores them
245: in a two element array.
247: Output Parameter:
248: + available - PETSC_TRUE if the val returned is valid
249: - val - the norm
251: Notes:
252: $ NORM_1 denotes sum_i |x_i|
253: $ NORM_2 denotes sqrt(sum_i (x_i)^2)
254: $ NORM_INFINITY denotes max_i |x_i|
256: Level: intermediate
258: Performance Issues:
259: $ per-processor memory bandwidth
260: $ interprocessor latency
261: $ work load inbalance that causes certain processes to arrive much earlier than others
263: Compile Option:
264: PETSC_HAVE_SLOW_BLAS_NORM2 will cause a C (loop unrolled) version of the norm to be used, rather
265: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
266: (as opposed to vendor provided) because the FORTRAN BLAS NRM2() routine is very slow.
268: Concepts: norm
269: Concepts: vector^norm
271: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), VecNorm()
272: VecNormBegin(), VecNormEnd()
274: @*/
275: PetscErrorCode VecNormAvailable(Vec x,NormType type,PetscBool *available,PetscReal *val)
276: {
284: *available = PETSC_FALSE;
285: if (type!=NORM_1_AND_2) {
286: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,*available);
287: }
288: return(0);
289: }
291: /*@
292: VecNormalize - Normalizes a vector by 2-norm.
294: Collective on Vec
296: Input Parameters:
297: + x - the vector
299: Output Parameter:
300: . x - the normalized vector
301: - val - the vector norm before normalization
303: Level: intermediate
305: Concepts: vector^normalizing
306: Concepts: normalizing^vector
308: @*/
309: PetscErrorCode VecNormalize(Vec x,PetscReal *val)
310: {
312: PetscReal norm;
317: PetscLogEventBegin(VEC_Normalize,x,0,0,0);
318: VecNorm(x,NORM_2,&norm);
319: if (norm == 0.0) {
320: PetscInfo(x,"Vector of zero norm can not be normalized; Returning only the zero norm\n");
321: } else if (norm != 1.0) {
322: PetscScalar tmp = 1.0/norm;
323: VecScale(x,tmp);
324: }
325: if (val) *val = norm;
326: PetscLogEventEnd(VEC_Normalize,x,0,0,0);
327: return(0);
328: }
330: /*@C
331: VecMax - Determines the vector component with maximum real part and its location.
333: Collective on Vec
335: Input Parameter:
336: . x - the vector
338: Output Parameters:
339: + p - the location of val (pass NULL if you don't want this)
340: - val - the maximum component
342: Notes:
343: Returns the value PETSC_MIN_REAL and p = -1 if the vector is of length 0.
345: Returns the smallest index with the maximum value
346: Level: intermediate
348: Concepts: maximum^of vector
349: Concepts: vector^maximum value
351: .seealso: VecNorm(), VecMin()
352: @*/
353: PetscErrorCode VecMax(Vec x,PetscInt *p,PetscReal *val)
354: {
361: PetscLogEventBegin(VEC_Max,x,0,0,0);
362: (*x->ops->max)(x,p,val);
363: PetscLogEventEnd(VEC_Max,x,0,0,0);
364: return(0);
365: }
367: /*@C
368: VecMin - Determines the vector component with minimum real part and its location.
370: Collective on Vec
372: Input Parameters:
373: . x - the vector
375: Output Parameter:
376: + p - the location of val (pass NULL if you don't want this location)
377: - val - the minimum component
379: Level: intermediate
381: Notes:
382: Returns the value PETSC_MAX_REAL and p = -1 if the vector is of length 0.
384: This returns the smallest index with the minumum value
386: Concepts: minimum^of vector
387: Concepts: vector^minimum entry
389: .seealso: VecMax()
390: @*/
391: PetscErrorCode VecMin(Vec x,PetscInt *p,PetscReal *val)
392: {
399: PetscLogEventBegin(VEC_Min,x,0,0,0);
400: (*x->ops->min)(x,p,val);
401: PetscLogEventEnd(VEC_Min,x,0,0,0);
402: return(0);
403: }
405: /*@
406: VecTDot - Computes an indefinite vector dot product. That is, this
407: routine does NOT use the complex conjugate.
409: Collective on Vec
411: Input Parameters:
412: . x, y - the vectors
414: Output Parameter:
415: . val - the dot product
417: Notes for Users of Complex Numbers:
418: For complex vectors, VecTDot() computes the indefinite form
419: $ val = (x,y) = y^T x,
420: where y^T denotes the transpose of y.
422: Use VecDot() for the inner product
423: $ val = (x,y) = y^H x,
424: where y^H denotes the conjugate transpose of y.
426: Level: intermediate
428: Concepts: inner product^non-Hermitian
429: Concepts: vector^inner product
430: Concepts: non-Hermitian inner product
432: .seealso: VecDot(), VecMTDot()
433: @*/
434: PetscErrorCode VecTDot(Vec x,Vec y,PetscScalar *val)
435: {
445: VecCheckSameSize(x,1,y,2);
447: PetscLogEventBegin(VEC_TDot,x,y,0,0);
448: (*x->ops->tdot)(x,y,val);
449: PetscLogEventEnd(VEC_TDot,x,y,0,0);
450: return(0);
451: }
453: /*@
454: VecScale - Scales a vector.
456: Not collective on Vec
458: Input Parameters:
459: + x - the vector
460: - alpha - the scalar
462: Output Parameter:
463: . x - the scaled vector
465: Note:
466: For a vector with n components, VecScale() computes
467: $ x[i] = alpha * x[i], for i=1,...,n.
469: Level: intermediate
471: Concepts: vector^scaling
472: Concepts: scaling^vector
474: @*/
475: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
476: {
477: PetscReal norms[4] = {0.0,0.0,0.0, 0.0};
478: PetscBool flgs[4];
480: PetscInt i;
485: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
486: PetscLogEventBegin(VEC_Scale,x,0,0,0);
487: if (alpha != (PetscScalar)1.0) {
488: /* get current stashed norms */
489: for (i=0; i<4; i++) {
490: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
491: }
492: (*x->ops->scale)(x,alpha);
493: PetscObjectStateIncrease((PetscObject)x);
494: /* put the scaled stashed norms back into the Vec */
495: for (i=0; i<4; i++) {
496: if (flgs[i]) {
497: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],PetscAbsScalar(alpha)*norms[i]);
498: }
499: }
500: }
501: PetscLogEventEnd(VEC_Scale,x,0,0,0);
502: return(0);
503: }
505: /*@
506: VecSet - Sets all components of a vector to a single scalar value.
508: Logically Collective on Vec
510: Input Parameters:
511: + x - the vector
512: - alpha - the scalar
514: Output Parameter:
515: . x - the vector
517: Note:
518: For a vector of dimension n, VecSet() computes
519: $ x[i] = alpha, for i=1,...,n,
520: so that all vector entries then equal the identical
521: scalar value, alpha. Use the more general routine
522: VecSetValues() to set different vector entries.
524: You CANNOT call this after you have called VecSetValues() but before you call
525: VecAssemblyBegin/End().
527: Level: beginner
529: .seealso VecSetValues(), VecSetValuesBlocked(), VecSetRandom()
531: Concepts: vector^setting to constant
533: @*/
534: PetscErrorCode VecSet(Vec x,PetscScalar alpha)
535: {
536: PetscReal val;
542: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"You cannot call this after you have called VecSetValues() but\n before you have called VecAssemblyBegin/End()");
544: VecSetErrorIfLocked(x,1);
546: PetscLogEventBegin(VEC_Set,x,0,0,0);
547: (*x->ops->set)(x,alpha);
548: PetscLogEventEnd(VEC_Set,x,0,0,0);
549: PetscObjectStateIncrease((PetscObject)x);
551: /* norms can be simply set (if |alpha|*N not too large) */
552: val = PetscAbsScalar(alpha);
553: if (x->map->N == 0) {
554: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],0.0l);
555: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],0.0);
556: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],0.0);
557: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],0.0);
558: } else if (val > PETSC_MAX_REAL/x->map->N) {
559: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
560: } else {
561: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],x->map->N * val);
562: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
563: val = PetscSqrtReal((PetscReal)x->map->N) * val;
564: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],val);
565: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],val);
566: }
567: return(0);
568: }
571: /*@
572: VecAXPY - Computes y = alpha x + y.
574: Logically Collective on Vec
576: Input Parameters:
577: + alpha - the scalar
578: - x, y - the vectors
580: Output Parameter:
581: . y - output vector
583: Level: intermediate
585: Notes:
586: x and y MUST be different vectors
588: Concepts: vector^BLAS
589: Concepts: BLAS
591: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY()
592: @*/
593: PetscErrorCode VecAXPY(Vec y,PetscScalar alpha,Vec x)
594: {
603: VecCheckSameSize(x,1,y,3);
604: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
606: VecSetErrorIfLocked(y,1);
608: VecLockReadPush(x);
609: PetscLogEventBegin(VEC_AXPY,x,y,0,0);
610: (*y->ops->axpy)(y,alpha,x);
611: PetscLogEventEnd(VEC_AXPY,x,y,0,0);
612: VecLockReadPop(x);
613: PetscObjectStateIncrease((PetscObject)y);
614: return(0);
615: }
617: /*@
618: VecAXPBY - Computes y = alpha x + beta y.
620: Logically Collective on Vec
622: Input Parameters:
623: + alpha,beta - the scalars
624: - x, y - the vectors
626: Output Parameter:
627: . y - output vector
629: Level: intermediate
631: Notes:
632: x and y MUST be different vectors
634: Concepts: BLAS
635: Concepts: vector^BLAS
637: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
638: @*/
639: PetscErrorCode VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
640: {
649: VecCheckSameSize(x,1,y,4);
650: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y cannot be the same vector");
654: PetscLogEventBegin(VEC_AXPY,x,y,0,0);
655: (*y->ops->axpby)(y,alpha,beta,x);
656: PetscLogEventEnd(VEC_AXPY,x,y,0,0);
657: PetscObjectStateIncrease((PetscObject)y);
658: return(0);
659: }
661: /*@
662: VecAXPBYPCZ - Computes z = alpha x + beta y + gamma z
664: Logically Collective on Vec
666: Input Parameters:
667: + alpha,beta, gamma - the scalars
668: - x, y, z - the vectors
670: Output Parameter:
671: . z - output vector
673: Level: intermediate
675: Notes:
676: x, y and z must be different vectors
678: Developer Note: alpha = 1 or gamma = 1 or gamma = 0.0 are handled as special cases
680: Concepts: BLAS
681: Concepts: vector^BLAS
683: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
684: @*/
685: PetscErrorCode VecAXPBYPCZ(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
686: {
698: VecCheckSameSize(x,1,y,5);
699: VecCheckSameSize(x,1,z,6);
700: if (x == y || x == z) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
701: if (y == z) SETERRQ(PetscObjectComm((PetscObject)y),PETSC_ERR_ARG_IDN,"x, y, and z must be different vectors");
706: PetscLogEventBegin(VEC_AXPBYPCZ,x,y,z,0);
707: (*y->ops->axpbypcz)(z,alpha,beta,gamma,x,y);
708: PetscLogEventEnd(VEC_AXPBYPCZ,x,y,z,0);
709: PetscObjectStateIncrease((PetscObject)z);
710: return(0);
711: }
713: /*@
714: VecAYPX - Computes y = x + alpha y.
716: Logically Collective on Vec
718: Input Parameters:
719: + alpha - the scalar
720: - x, y - the vectors
722: Output Parameter:
723: . y - output vector
725: Level: intermediate
727: Notes:
728: x and y MUST be different vectors
730: Concepts: vector^BLAS
731: Concepts: BLAS
733: .seealso: VecAXPY(), VecWAXPY()
734: @*/
735: PetscErrorCode VecAYPX(Vec y,PetscScalar alpha,Vec x)
736: {
745: VecCheckSameSize(x,1,y,3);
746: if (x == y) SETERRQ(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
749: PetscLogEventBegin(VEC_AYPX,x,y,0,0);
750: (*y->ops->aypx)(y,alpha,x);
751: PetscLogEventEnd(VEC_AYPX,x,y,0,0);
752: PetscObjectStateIncrease((PetscObject)y);
753: return(0);
754: }
757: /*@
758: VecWAXPY - Computes w = alpha x + y.
760: Logically Collective on Vec
762: Input Parameters:
763: + alpha - the scalar
764: - x, y - the vectors
766: Output Parameter:
767: . w - the result
769: Level: intermediate
771: Notes:
772: w cannot be either x or y, but x and y can be the same
774: Concepts: vector^BLAS
775: Concepts: BLAS
777: .seealso: VecAXPY(), VecAYPX(), VecAXPBY()
778: @*/
779: PetscErrorCode VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
780: {
792: VecCheckSameSize(x,3,y,4);
793: VecCheckSameSize(x,3,w,1);
794: if (w == y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector y, suggest VecAXPY()");
795: if (w == x) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Result vector w cannot be same as input vector x, suggest VecAYPX()");
798: PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
799: (*w->ops->waxpy)(w,alpha,x,y);
800: PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
801: PetscObjectStateIncrease((PetscObject)w);
802: return(0);
803: }
806: /*@C
807: VecSetValues - Inserts or adds values into certain locations of a vector.
809: Not Collective
811: Input Parameters:
812: + x - vector to insert in
813: . ni - number of elements to add
814: . ix - indices where to add
815: . y - array of values
816: - iora - either INSERT_VALUES or ADD_VALUES, where
817: ADD_VALUES adds values to any existing entries, and
818: INSERT_VALUES replaces existing entries with new values
820: Notes:
821: VecSetValues() sets x[ix[i]] = y[i], for i=0,...,ni-1.
823: Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
824: options cannot be mixed without intervening calls to the assembly
825: routines.
827: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
828: MUST be called after all calls to VecSetValues() have been completed.
830: VecSetValues() uses 0-based indices in Fortran as well as in C.
832: If you call VecSetOption(x, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE),
833: negative indices may be passed in ix. These rows are
834: simply ignored. This allows easily inserting element load matrices
835: with homogeneous Dirchlet boundary conditions that you don't want represented
836: in the vector.
838: Level: beginner
840: Concepts: vector^setting values
842: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
843: VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
844: @*/
845: PetscErrorCode VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
846: {
851: if (!ni) return(0);
855: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
856: (*x->ops->setvalues)(x,ni,ix,y,iora);
857: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
858: PetscObjectStateIncrease((PetscObject)x);
859: return(0);
860: }
862: /*@
863: VecGetValues - Gets values from certain locations of a vector. Currently
864: can only get values on the same processor
866: Not Collective
868: Input Parameters:
869: + x - vector to get values from
870: . ni - number of elements to get
871: - ix - indices where to get them from (in global 1d numbering)
873: Output Parameter:
874: . y - array of values
876: Notes:
877: The user provides the allocated array y; it is NOT allocated in this routine
879: VecGetValues() gets y[i] = x[ix[i]], for i=0,...,ni-1.
881: VecAssemblyBegin() and VecAssemblyEnd() MUST be called before calling this
883: VecGetValues() uses 0-based indices in Fortran as well as in C.
885: If you call VecSetOption(x, VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE),
886: negative indices may be passed in ix. These rows are
887: simply ignored.
889: Level: beginner
891: Concepts: vector^getting values
893: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues()
894: @*/
895: PetscErrorCode VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
896: {
901: if (!ni) return(0);
905: (*x->ops->getvalues)(x,ni,ix,y);
906: return(0);
907: }
909: /*@C
910: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
912: Not Collective
914: Input Parameters:
915: + x - vector to insert in
916: . ni - number of blocks to add
917: . ix - indices where to add in block count, rather than element count
918: . y - array of values
919: - iora - either INSERT_VALUES or ADD_VALUES, where
920: ADD_VALUES adds values to any existing entries, and
921: INSERT_VALUES replaces existing entries with new values
923: Notes:
924: VecSetValuesBlocked() sets x[bs*ix[i]+j] = y[bs*i+j],
925: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
927: Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
928: options cannot be mixed without intervening calls to the assembly
929: routines.
931: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
932: MUST be called after all calls to VecSetValuesBlocked() have been completed.
934: VecSetValuesBlocked() uses 0-based indices in Fortran as well as in C.
936: Negative indices may be passed in ix, these rows are
937: simply ignored. This allows easily inserting element load matrices
938: with homogeneous Dirchlet boundary conditions that you don't want represented
939: in the vector.
941: Level: intermediate
943: Concepts: vector^setting values blocked
945: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
946: VecSetValues()
947: @*/
948: PetscErrorCode VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
949: {
954: if (!ni) return(0);
958: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
959: (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
960: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
961: PetscObjectStateIncrease((PetscObject)x);
962: return(0);
963: }
966: /*@C
967: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
968: using a local ordering of the nodes.
970: Not Collective
972: Input Parameters:
973: + x - vector to insert in
974: . ni - number of elements to add
975: . ix - indices where to add
976: . y - array of values
977: - iora - either INSERT_VALUES or ADD_VALUES, where
978: ADD_VALUES adds values to any existing entries, and
979: INSERT_VALUES replaces existing entries with new values
981: Level: intermediate
983: Notes:
984: VecSetValuesLocal() sets x[ix[i]] = y[i], for i=0,...,ni-1.
986: Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
987: options cannot be mixed without intervening calls to the assembly
988: routines.
990: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
991: MUST be called after all calls to VecSetValuesLocal() have been completed.
993: VecSetValuesLocal() uses 0-based indices in Fortran as well as in C.
995: Concepts: vector^setting values with local numbering
997: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetLocalToGlobalMapping(),
998: VecSetValuesBlockedLocal()
999: @*/
1000: PetscErrorCode VecSetValuesLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1001: {
1003: PetscInt lixp[128],*lix = lixp;
1007: if (!ni) return(0);
1012: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1013: if (!x->ops->setvalueslocal) {
1014: if (!x->map->mapping) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMapping()");
1015: if (ni > 128) {
1016: PetscMalloc1(ni,&lix);
1017: }
1018: ISLocalToGlobalMappingApply(x->map->mapping,ni,(PetscInt*)ix,lix);
1019: (*x->ops->setvalues)(x,ni,lix,y,iora);
1020: if (ni > 128) {
1021: PetscFree(lix);
1022: }
1023: } else {
1024: (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
1025: }
1026: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1027: PetscObjectStateIncrease((PetscObject)x);
1028: return(0);
1029: }
1031: /*@
1032: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1033: using a local ordering of the nodes.
1035: Not Collective
1037: Input Parameters:
1038: + x - vector to insert in
1039: . ni - number of blocks to add
1040: . ix - indices where to add in block count, not element count
1041: . y - array of values
1042: - iora - either INSERT_VALUES or ADD_VALUES, where
1043: ADD_VALUES adds values to any existing entries, and
1044: INSERT_VALUES replaces existing entries with new values
1046: Level: intermediate
1048: Notes:
1049: VecSetValuesBlockedLocal() sets x[bs*ix[i]+j] = y[bs*i+j],
1050: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with VecSetBlockSize().
1052: Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
1053: options cannot be mixed without intervening calls to the assembly
1054: routines.
1056: These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
1057: MUST be called after all calls to VecSetValuesBlockedLocal() have been completed.
1059: VecSetValuesBlockedLocal() uses 0-based indices in Fortran as well as in C.
1062: Concepts: vector^setting values blocked with local numbering
1064: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
1065: VecSetLocalToGlobalMapping()
1066: @*/
1067: PetscErrorCode VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
1068: {
1070: PetscInt lixp[128],*lix = lixp;
1074: if (!ni) return(0);
1078: if (ni > 128) {
1079: PetscMalloc1(ni,&lix);
1080: }
1082: PetscLogEventBegin(VEC_SetValues,x,0,0,0);
1083: ISLocalToGlobalMappingApplyBlock(x->map->mapping,ni,(PetscInt*)ix,lix);
1084: (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
1085: PetscLogEventEnd(VEC_SetValues,x,0,0,0);
1086: if (ni > 128) {
1087: PetscFree(lix);
1088: }
1089: PetscObjectStateIncrease((PetscObject)x);
1090: return(0);
1091: }
1093: /*@
1094: VecMTDot - Computes indefinite vector multiple dot products.
1095: That is, it does NOT use the complex conjugate.
1097: Collective on Vec
1099: Input Parameters:
1100: + x - one vector
1101: . nv - number of vectors
1102: - y - array of vectors. Note that vectors are pointers
1104: Output Parameter:
1105: . val - array of the dot products
1107: Notes for Users of Complex Numbers:
1108: For complex vectors, VecMTDot() computes the indefinite form
1109: $ val = (x,y) = y^T x,
1110: where y^T denotes the transpose of y.
1112: Use VecMDot() for the inner product
1113: $ val = (x,y) = y^H x,
1114: where y^H denotes the conjugate transpose of y.
1116: Level: intermediate
1118: Concepts: inner product^multiple
1119: Concepts: vector^multiple inner products
1121: .seealso: VecMDot(), VecTDot()
1122: @*/
1123: PetscErrorCode VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1124: {
1130: if (!nv) return(0);
1137: VecCheckSameSize(x,1,*y,3);
1139: PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
1140: (*x->ops->mtdot)(x,nv,y,val);
1141: PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
1142: return(0);
1143: }
1145: /*@
1146: VecMDot - Computes vector multiple dot products.
1148: Collective on Vec
1150: Input Parameters:
1151: + x - one vector
1152: . nv - number of vectors
1153: - y - array of vectors.
1155: Output Parameter:
1156: . val - array of the dot products (does not allocate the array)
1158: Notes for Users of Complex Numbers:
1159: For complex vectors, VecMDot() computes
1160: $ val = (x,y) = y^H x,
1161: where y^H denotes the conjugate transpose of y.
1163: Use VecMTDot() for the indefinite form
1164: $ val = (x,y) = y^T x,
1165: where y^T denotes the transpose of y.
1167: Level: intermediate
1169: Concepts: inner product^multiple
1170: Concepts: vector^multiple inner products
1172: .seealso: VecMTDot(), VecDot()
1173: @*/
1174: PetscErrorCode VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar val[])
1175: {
1181: if (!nv) return(0);
1182: if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1189: VecCheckSameSize(x,1,*y,3);
1191: PetscLogEventBegin(VEC_MDot,x,*y,0,0);
1192: (*x->ops->mdot)(x,nv,y,val);
1193: PetscLogEventEnd(VEC_MDot,x,*y,0,0);
1194: return(0);
1195: }
1197: /*@
1198: VecMAXPY - Computes y = y + sum alpha[j] x[j]
1200: Logically Collective on Vec
1202: Input Parameters:
1203: + nv - number of scalars and x-vectors
1204: . alpha - array of scalars
1205: . y - one vector
1206: - x - array of vectors
1208: Level: intermediate
1210: Notes:
1211: y cannot be any of the x vectors
1213: Concepts: BLAS
1215: .seealso: VecAXPY(), VecWAXPY(), VecAYPX()
1216: @*/
1217: PetscErrorCode VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec x[])
1218: {
1220: PetscInt i;
1225: if (!nv) return(0);
1226: if (nv < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",nv);
1233: VecCheckSameSize(y,1,*x,4);
1236: PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
1237: (*y->ops->maxpy)(y,nv,alpha,x);
1238: PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
1239: PetscObjectStateIncrease((PetscObject)y);
1240: return(0);
1241: }
1243: /*@
1244: VecGetSubVector - Gets a vector representing part of another vector
1246: Collective on IS (and Vec if nonlocal entries are needed)
1248: Input Arguments:
1249: + X - vector from which to extract a subvector
1250: - is - index set representing portion of X to extract
1252: Output Arguments:
1253: . Y - subvector corresponding to is
1255: Level: advanced
1257: Notes:
1258: The subvector Y should be returned with VecRestoreSubVector().
1260: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1261: modifying the subvector. Other non-overlapping subvectors can still be obtained from X using this function.
1263: .seealso: MatCreateSubMatrix()
1264: @*/
1265: PetscErrorCode VecGetSubVector(Vec X,IS is,Vec *Y)
1266: {
1267: PetscErrorCode ierr;
1268: Vec Z;
1274: if (X->ops->getsubvector) {
1275: (*X->ops->getsubvector)(X,is,&Z);
1276: } else { /* Default implementation currently does no caching */
1277: PetscInt gstart,gend,start;
1278: PetscBool contiguous,gcontiguous;
1279: VecGetOwnershipRange(X,&gstart,&gend);
1280: ISContiguousLocal(is,gstart,gend,&start,&contiguous);
1281: MPIU_Allreduce(&contiguous,&gcontiguous,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)is));
1282: if (gcontiguous) { /* We can do a no-copy implementation */
1283: PetscInt n,N,bs;
1284: PetscInt state;
1286: ISGetSize(is,&N);
1287: ISGetLocalSize(is,&n);
1288: VecGetBlockSize(X,&bs);
1289: if (n%bs || bs == 1 || !n) bs = -1; /* Do not decide block size if we do not have to */
1290: VecLockGet(X,&state);
1291: if (state) {
1292: const PetscScalar *x;
1293: VecGetArrayRead(X,&x);
1294: VecCreate(PetscObjectComm((PetscObject)X),&Z);
1295: VecSetType(Z,((PetscObject)X)->type_name);
1296: VecSetSizes(Z,n,N);
1297: VecSetBlockSize(Z,bs);
1298: VecPlaceArray(Z,(PetscScalar*)x+start);
1299: VecLockReadPush(Z);
1300: VecRestoreArrayRead(X,&x);
1301: } else {
1302: PetscScalar *x;
1303: VecGetArray(X,&x);
1304: VecCreate(PetscObjectComm((PetscObject)X),&Z);
1305: VecSetType(Z,((PetscObject)X)->type_name);
1306: VecSetSizes(Z,n,N);
1307: VecSetBlockSize(Z,bs);
1308: VecPlaceArray(Z,(PetscScalar*)x+start);
1309: VecRestoreArray(X,&x);
1310: }
1311: } else { /* Have to create a scatter and do a copy */
1312: VecScatter scatter;
1313: PetscInt n,N;
1314: ISGetLocalSize(is,&n);
1315: ISGetSize(is,&N);
1316: VecCreate(PetscObjectComm((PetscObject)is),&Z);
1317: VecSetSizes(Z,n,N);
1318: VecSetType(Z,((PetscObject)X)->type_name);
1319: VecScatterCreate(X,is,Z,NULL,&scatter);
1320: VecScatterBegin(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1321: VecScatterEnd(scatter,X,Z,INSERT_VALUES,SCATTER_FORWARD);
1322: PetscObjectCompose((PetscObject)Z,"VecGetSubVector_Scatter",(PetscObject)scatter);
1323: VecScatterDestroy(&scatter);
1324: }
1325: }
1326: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1327: if (VecGetSubVectorSavedStateId < 0) {PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId);}
1328: PetscObjectComposedDataSetInt((PetscObject)Z,VecGetSubVectorSavedStateId,1);
1329: *Y = Z;
1330: return(0);
1331: }
1333: /*@
1334: VecRestoreSubVector - Restores a subvector extracted using VecGetSubVector()
1336: Collective on IS (and Vec if nonlocal entries need to be written)
1338: Input Arguments:
1339: + X - vector from which subvector was obtained
1340: . is - index set representing the subset of X
1341: - Y - subvector being restored
1343: Level: advanced
1345: .seealso: VecGetSubVector()
1346: @*/
1347: PetscErrorCode VecRestoreSubVector(Vec X,IS is,Vec *Y)
1348: {
1356: if (X->ops->restoresubvector) {
1357: (*X->ops->restoresubvector)(X,is,Y);
1358: } else {
1359: PETSC_UNUSED PetscObjectState dummystate = 0;
1360: PetscBool valid;
1361: PetscObjectComposedDataGetInt((PetscObject)*Y,VecGetSubVectorSavedStateId,dummystate,valid);
1362: if (!valid) {
1363: VecScatter scatter;
1365: PetscObjectQuery((PetscObject)*Y,"VecGetSubVector_Scatter",(PetscObject*)&scatter);
1366: if (scatter) {
1367: VecScatterBegin(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1368: VecScatterEnd(scatter,*Y,X,INSERT_VALUES,SCATTER_REVERSE);
1369: }
1370: }
1371: VecDestroy(Y);
1372: }
1373: return(0);
1374: }
1376: /*@
1377: VecGetLocalVectorRead - Maps the local portion of a vector into a
1378: vector. You must call VecRestoreLocalVectorRead() when the local
1379: vector is no longer needed.
1381: Not collective.
1383: Input parameter:
1384: . v - The vector for which the local vector is desired.
1386: Output parameter:
1387: . w - Upon exit this contains the local vector.
1389: Level: beginner
1391: Notes:
1392: This function is similar to VecGetArrayRead() which maps the local
1393: portion into a raw pointer. VecGetLocalVectorRead() is usually
1394: almost as efficient as VecGetArrayRead() but in certain circumstances
1395: VecGetLocalVectorRead() can be much more efficient than
1396: VecGetArrayRead(). This is because the construction of a contiguous
1397: array representing the vector data required by VecGetArrayRead() can
1398: be an expensive operation for certain vector types. For example, for
1399: GPU vectors VecGetArrayRead() requires that the data between device
1400: and host is synchronized.
1402: Unlike VecGetLocalVector(), this routine is not collective and
1403: preserves cached information.
1405: .seealso: VecRestoreLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1406: @*/
1407: PetscErrorCode VecGetLocalVectorRead(Vec v,Vec w)
1408: {
1410: PetscScalar *a;
1415: VecCheckSameLocalSize(v,1,w,2);
1416: if (v->ops->getlocalvectorread) {
1417: (*v->ops->getlocalvectorread)(v,w);
1418: } else {
1419: VecGetArrayRead(v,(const PetscScalar**)&a);
1420: VecPlaceArray(w,a);
1421: }
1422: return(0);
1423: }
1425: /*@
1426: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1427: previously mapped into a vector using VecGetLocalVectorRead().
1429: Not collective.
1431: Input parameter:
1432: . v - The local portion of this vector was previously mapped into w using VecGetLocalVectorRead().
1433: . w - The vector into which the local portion of v was mapped.
1435: Level: beginner
1437: .seealso: VecGetLocalVectorRead(), VecGetLocalVector(), VecGetArrayRead(), VecGetArray()
1438: @*/
1439: PetscErrorCode VecRestoreLocalVectorRead(Vec v,Vec w)
1440: {
1442: PetscScalar *a;
1447: if (v->ops->restorelocalvectorread) {
1448: (*v->ops->restorelocalvectorread)(v,w);
1449: } else {
1450: VecGetArrayRead(w,(const PetscScalar**)&a);
1451: VecRestoreArrayRead(v,(const PetscScalar**)&a);
1452: VecResetArray(w);
1453: }
1454: return(0);
1455: }
1457: /*@
1458: VecGetLocalVector - Maps the local portion of a vector into a
1459: vector.
1461: Collective on v, not collective on w.
1463: Input parameter:
1464: . v - The vector for which the local vector is desired.
1466: Output parameter:
1467: . w - Upon exit this contains the local vector.
1469: Level: beginner
1471: Notes:
1472: This function is similar to VecGetArray() which maps the local
1473: portion into a raw pointer. VecGetLocalVector() is usually about as
1474: efficient as VecGetArray() but in certain circumstances
1475: VecGetLocalVector() can be much more efficient than VecGetArray().
1476: This is because the construction of a contiguous array representing
1477: the vector data required by VecGetArray() can be an expensive
1478: operation for certain vector types. For example, for GPU vectors
1479: VecGetArray() requires that the data between device and host is
1480: synchronized.
1482: .seealso: VecRestoreLocalVector(), VecGetLocalVectorRead(), VecGetArrayRead(), VecGetArray()
1483: @*/
1484: PetscErrorCode VecGetLocalVector(Vec v,Vec w)
1485: {
1487: PetscScalar *a;
1492: VecCheckSameLocalSize(v,1,w,2);
1493: if (v->ops->getlocalvector) {
1494: (*v->ops->getlocalvector)(v,w);
1495: } else {
1496: VecGetArray(v,&a);
1497: VecPlaceArray(w,a);
1498: }
1499: return(0);
1500: }
1502: /*@
1503: VecRestoreLocalVector - Unmaps the local portion of a vector
1504: previously mapped into a vector using VecGetLocalVector().
1506: Logically collective.
1508: Input parameter:
1509: . v - The local portion of this vector was previously mapped into w using VecGetLocalVector().
1510: . w - The vector into which the local portion of v was mapped.
1512: Level: beginner
1514: .seealso: VecGetLocalVector(), VecGetLocalVectorRead(), VecRestoreLocalVectorRead(), LocalVectorRead(), VecGetArrayRead(), VecGetArray()
1515: @*/
1516: PetscErrorCode VecRestoreLocalVector(Vec v,Vec w)
1517: {
1519: PetscScalar *a;
1524: if (v->ops->restorelocalvector) {
1525: (*v->ops->restorelocalvector)(v,w);
1526: } else {
1527: VecGetArray(w,&a);
1528: VecRestoreArray(v,&a);
1529: VecResetArray(w);
1530: }
1531: return(0);
1532: }
1534: /*@C
1535: VecGetArray - Returns a pointer to a contiguous array that contains this
1536: processor's portion of the vector data. For the standard PETSc
1537: vectors, VecGetArray() returns a pointer to the local data array and
1538: does not use any copies. If the underlying vector data is not stored
1539: in a contiguous array this routine will copy the data to a contiguous
1540: array and return a pointer to that. You MUST call VecRestoreArray()
1541: when you no longer need access to the array.
1543: Logically Collective on Vec
1545: Input Parameter:
1546: . x - the vector
1548: Output Parameter:
1549: . a - location to put pointer to the array
1551: Fortran Note:
1552: This routine is used differently from Fortran 77
1553: $ Vec x
1554: $ PetscScalar x_array(1)
1555: $ PetscOffset i_x
1556: $ PetscErrorCode ierr
1557: $ call VecGetArray(x,x_array,i_x,ierr)
1558: $
1559: $ Access first local entry in vector with
1560: $ value = x_array(i_x + 1)
1561: $
1562: $ ...... other code
1563: $ call VecRestoreArray(x,x_array,i_x,ierr)
1564: For Fortran 90 see VecGetArrayF90()
1566: See the Fortran chapter of the users manual and
1567: petsc/src/snes/examples/tutorials/ex5f.F for details.
1569: Level: beginner
1571: Concepts: vector^accessing local values
1573: .seealso: VecRestoreArray(), VecGetArrayRead(), VecGetArrays(), VecGetArrayF90(), VecGetArrayReadF90(), VecPlaceArray(), VecGetArray2d(),
1574: VecGetArrayPair(), VecRestoreArrayPair()
1575: @*/
1576: PetscErrorCode VecGetArray(Vec x,PetscScalar **a)
1577: {
1579: #if defined(PETSC_HAVE_VIENNACL)
1580: PetscBool is_viennacltype = PETSC_FALSE;
1581: #endif
1585: VecSetErrorIfLocked(x,1);
1586: if (x->petscnative) {
1587: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1588: if (x->valid_GPU_array == PETSC_OFFLOAD_GPU) {
1589: #if defined(PETSC_HAVE_VIENNACL)
1590: PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1591: if (is_viennacltype) {
1592: VecViennaCLCopyFromGPU(x);
1593: } else
1594: #endif
1595: {
1596: #if defined(PETSC_HAVE_CUDA)
1597: VecCUDACopyFromGPU(x);
1598: #endif
1599: }
1600: } else if (x->valid_GPU_array == PETSC_OFFLOAD_UNALLOCATED) {
1601: #if defined(PETSC_HAVE_VIENNACL)
1602: PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1603: if (is_viennacltype) {
1604: VecViennaCLAllocateCheckHost(x);
1605: } else
1606: #endif
1607: {
1608: #if defined(PETSC_HAVE_CUDA)
1609: VecCUDAAllocateCheckHost(x);
1610: #endif
1611: }
1612: }
1613: #endif
1614: *a = *((PetscScalar**)x->data);
1615: } else {
1616: if (x->ops->getarray) {
1617: (*x->ops->getarray)(x,a);
1618: } else SETERRQ1(PetscObjectComm((PetscObject)x),PETSC_ERR_SUP,"Cannot get array for vector type \"%s\"",((PetscObject)x)->type_name);
1619: }
1620: return(0);
1621: }
1623: /*@C
1624: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
1626: Not Collective
1628: Input Parameters:
1629: . x - the vector
1631: Output Parameter:
1632: . a - the array
1634: Level: beginner
1636: Notes:
1637: The array must be returned using a matching call to VecRestoreArrayRead().
1639: Unlike VecGetArray(), this routine is not collective and preserves cached information like vector norms.
1641: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
1642: implementations may require a copy, but must such implementations should cache the contiguous representation so that
1643: only one copy is performed when this routine is called multiple times in sequence.
1645: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1646: @*/
1647: PetscErrorCode VecGetArrayRead(Vec x,const PetscScalar **a)
1648: {
1650: #if defined(PETSC_HAVE_VIENNACL)
1651: PetscBool is_viennacltype = PETSC_FALSE;
1652: #endif
1656: if (x->petscnative) {
1657: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1658: if (x->valid_GPU_array == PETSC_OFFLOAD_GPU) {
1659: #if defined(PETSC_HAVE_VIENNACL)
1660: PetscObjectTypeCompareAny((PetscObject)x,&is_viennacltype,VECSEQVIENNACL,VECMPIVIENNACL,VECVIENNACL,"");
1661: if (is_viennacltype) {
1662: VecViennaCLCopyFromGPU(x);
1663: } else
1664: #endif
1665: {
1666: #if defined(PETSC_HAVE_CUDA)
1667: VecCUDACopyFromGPU(x);
1668: #endif
1669: }
1670: }
1671: #endif
1672: *a = *((PetscScalar **)x->data);
1673: } else if (x->ops->getarrayread) {
1674: (*x->ops->getarrayread)(x,a);
1675: } else {
1676: (*x->ops->getarray)(x,(PetscScalar**)a);
1677: }
1678: return(0);
1679: }
1681: /*@C
1682: VecGetArrays - Returns a pointer to the arrays in a set of vectors
1683: that were created by a call to VecDuplicateVecs(). You MUST call
1684: VecRestoreArrays() when you no longer need access to the array.
1686: Logically Collective on Vec
1688: Input Parameter:
1689: + x - the vectors
1690: - n - the number of vectors
1692: Output Parameter:
1693: . a - location to put pointer to the array
1695: Fortran Note:
1696: This routine is not supported in Fortran.
1698: Level: intermediate
1700: .seealso: VecGetArray(), VecRestoreArrays()
1701: @*/
1702: PetscErrorCode VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1703: {
1705: PetscInt i;
1706: PetscScalar **q;
1712: if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1713: PetscMalloc1(n,&q);
1714: for (i=0; i<n; ++i) {
1715: VecGetArray(x[i],&q[i]);
1716: }
1717: *a = q;
1718: return(0);
1719: }
1721: /*@C
1722: VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1723: has been called.
1725: Logically Collective on Vec
1727: Input Parameters:
1728: + x - the vector
1729: . n - the number of vectors
1730: - a - location of pointer to arrays obtained from VecGetArrays()
1732: Notes:
1733: For regular PETSc vectors this routine does not involve any copies. For
1734: any special vectors that do not store local vector data in a contiguous
1735: array, this routine will copy the data back into the underlying
1736: vector data structure from the arrays obtained with VecGetArrays().
1738: Fortran Note:
1739: This routine is not supported in Fortran.
1741: Level: intermediate
1743: .seealso: VecGetArrays(), VecRestoreArray()
1744: @*/
1745: PetscErrorCode VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1746: {
1748: PetscInt i;
1749: PetscScalar **q = *a;
1756: for (i=0; i<n; ++i) {
1757: VecRestoreArray(x[i],&q[i]);
1758: }
1759: PetscFree(q);
1760: return(0);
1761: }
1763: /*@C
1764: VecRestoreArray - Restores a vector after VecGetArray() has been called.
1766: Logically Collective on Vec
1768: Input Parameters:
1769: + x - the vector
1770: - a - location of pointer to array obtained from VecGetArray()
1772: Level: beginner
1774: Notes:
1775: For regular PETSc vectors this routine does not involve any copies. For
1776: any special vectors that do not store local vector data in a contiguous
1777: array, this routine will copy the data back into the underlying
1778: vector data structure from the array obtained with VecGetArray().
1780: This routine actually zeros out the a pointer. This is to prevent accidental
1781: use of the array after it has been restored. If you pass null for a it will
1782: not zero the array pointer a.
1784: Fortran Note:
1785: This routine is used differently from Fortran 77
1786: $ Vec x
1787: $ PetscScalar x_array(1)
1788: $ PetscOffset i_x
1789: $ PetscErrorCode ierr
1790: $ call VecGetArray(x,x_array,i_x,ierr)
1791: $
1792: $ Access first local entry in vector with
1793: $ value = x_array(i_x + 1)
1794: $
1795: $ ...... other code
1796: $ call VecRestoreArray(x,x_array,i_x,ierr)
1798: See the Fortran chapter of the users manual and
1799: petsc/src/snes/examples/tutorials/ex5f.F for details.
1800: For Fortran 90 see VecRestoreArrayF90()
1802: .seealso: VecGetArray(), VecRestoreArrayRead(), VecRestoreArrays(), VecRestoreArrayF90(), VecRestoreArrayReadF90(), VecPlaceArray(), VecRestoreArray2d(),
1803: VecGetArrayPair(), VecRestoreArrayPair()
1804: @*/
1805: PetscErrorCode VecRestoreArray(Vec x,PetscScalar **a)
1806: {
1811: if (x->petscnative) {
1812: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1813: x->valid_GPU_array = PETSC_OFFLOAD_CPU;
1814: #endif
1815: } else {
1816: (*x->ops->restorearray)(x,a);
1817: }
1818: if (a) *a = NULL;
1819: PetscObjectStateIncrease((PetscObject)x);
1820: return(0);
1821: }
1823: /*@C
1824: VecRestoreArrayRead - Restore array obtained with VecGetArrayRead()
1826: Not Collective
1828: Input Parameters:
1829: + vec - the vector
1830: - array - the array
1832: Level: beginner
1834: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrayPair(), VecRestoreArrayPair()
1835: @*/
1836: PetscErrorCode VecRestoreArrayRead(Vec x,const PetscScalar **a)
1837: {
1842: if (x->petscnative) {
1843: /* nothing */
1844: } else if (x->ops->restorearrayread) {
1845: (*x->ops->restorearrayread)(x,a);
1846: } else {
1847: (*x->ops->restorearray)(x,(PetscScalar**)a);
1848: }
1849: if (a) *a = NULL;
1850: return(0);
1851: }
1853: /*@
1854: VecPlaceArray - Allows one to replace the array in a vector with an
1855: array provided by the user. This is useful to avoid copying an array
1856: into a vector.
1858: Not Collective
1860: Input Parameters:
1861: + vec - the vector
1862: - array - the array
1864: Notes:
1865: You can return to the original array with a call to VecResetArray()
1867: Level: developer
1869: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecResetArray()
1871: @*/
1872: PetscErrorCode VecPlaceArray(Vec vec,const PetscScalar array[])
1873: {
1880: if (vec->ops->placearray) {
1881: (*vec->ops->placearray)(vec,array);
1882: } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot place array in this type of vector");
1883: PetscObjectStateIncrease((PetscObject)vec);
1884: return(0);
1885: }
1887: /*@C
1888: VecReplaceArray - Allows one to replace the array in a vector with an
1889: array provided by the user. This is useful to avoid copying an array
1890: into a vector.
1892: Not Collective
1894: Input Parameters:
1895: + vec - the vector
1896: - array - the array
1898: Notes:
1899: This permanently replaces the array and frees the memory associated
1900: with the old array.
1902: The memory passed in MUST be obtained with PetscMalloc() and CANNOT be
1903: freed by the user. It will be freed when the vector is destroy.
1905: Not supported from Fortran
1907: Level: developer
1909: .seealso: VecGetArray(), VecRestoreArray(), VecPlaceArray(), VecResetArray()
1911: @*/
1912: PetscErrorCode VecReplaceArray(Vec vec,const PetscScalar array[])
1913: {
1919: if (vec->ops->replacearray) {
1920: (*vec->ops->replacearray)(vec,array);
1921: } else SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_SUP,"Cannot replace array in this type of vector");
1922: PetscObjectStateIncrease((PetscObject)vec);
1923: return(0);
1924: }
1926: /*MC
1927: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
1928: and makes them accessible via a Fortran90 pointer.
1930: Synopsis:
1931: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
1933: Collective on Vec
1935: Input Parameters:
1936: + x - a vector to mimic
1937: - n - the number of vectors to obtain
1939: Output Parameters:
1940: + y - Fortran90 pointer to the array of vectors
1941: - ierr - error code
1943: Example of Usage:
1944: .vb
1945: #include <petsc/finclude/petscvec.h>
1946: use petscvec
1948: Vec x
1949: Vec, pointer :: y(:)
1950: ....
1951: call VecDuplicateVecsF90(x,2,y,ierr)
1952: call VecSet(y(2),alpha,ierr)
1953: call VecSet(y(2),alpha,ierr)
1954: ....
1955: call VecDestroyVecsF90(2,y,ierr)
1956: .ve
1958: Notes:
1959: Not yet supported for all F90 compilers
1961: Use VecDestroyVecsF90() to free the space.
1963: Level: beginner
1965: .seealso: VecDestroyVecsF90(), VecDuplicateVecs()
1967: M*/
1969: /*MC
1970: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
1971: VecGetArrayF90().
1973: Synopsis:
1974: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
1976: Logically Collective on Vec
1978: Input Parameters:
1979: + x - vector
1980: - xx_v - the Fortran90 pointer to the array
1982: Output Parameter:
1983: . ierr - error code
1985: Example of Usage:
1986: .vb
1987: #include <petsc/finclude/petscvec.h>
1988: use petscvec
1990: PetscScalar, pointer :: xx_v(:)
1991: ....
1992: call VecGetArrayF90(x,xx_v,ierr)
1993: xx_v(3) = a
1994: call VecRestoreArrayF90(x,xx_v,ierr)
1995: .ve
1997: Level: beginner
1999: .seealso: VecGetArrayF90(), VecGetArray(), VecRestoreArray(), UsingFortran, VecRestoreArrayReadF90()
2001: M*/
2003: /*MC
2004: VecDestroyVecsF90 - Frees a block of vectors obtained with VecDuplicateVecsF90().
2006: Synopsis:
2007: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
2009: Collective on Vec
2011: Input Parameters:
2012: + n - the number of vectors previously obtained
2013: - x - pointer to array of vector pointers
2015: Output Parameter:
2016: . ierr - error code
2018: Notes:
2019: Not yet supported for all F90 compilers
2021: Level: beginner
2023: .seealso: VecDestroyVecs(), VecDuplicateVecsF90()
2025: M*/
2027: /*MC
2028: VecGetArrayF90 - Accesses a vector array from Fortran90. For default PETSc
2029: vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
2030: this routine is implementation dependent. You MUST call VecRestoreArrayF90()
2031: when you no longer need access to the array.
2033: Synopsis:
2034: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2036: Logically Collective on Vec
2038: Input Parameter:
2039: . x - vector
2041: Output Parameters:
2042: + xx_v - the Fortran90 pointer to the array
2043: - ierr - error code
2045: Example of Usage:
2046: .vb
2047: #include <petsc/finclude/petscvec.h>
2048: use petscvec
2050: PetscScalar, pointer :: xx_v(:)
2051: ....
2052: call VecGetArrayF90(x,xx_v,ierr)
2053: xx_v(3) = a
2054: call VecRestoreArrayF90(x,xx_v,ierr)
2055: .ve
2057: If you ONLY intend to read entries from the array and not change any entries you should use VecGetArrayReadF90().
2059: Level: beginner
2061: .seealso: VecRestoreArrayF90(), VecGetArray(), VecRestoreArray(), VecGetArrayReadF90(), UsingFortran
2063: M*/
2065: /*MC
2066: VecGetArrayReadF90 - Accesses a read only array from Fortran90. For default PETSc
2067: vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
2068: this routine is implementation dependent. You MUST call VecRestoreArrayReadF90()
2069: when you no longer need access to the array.
2071: Synopsis:
2072: VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2074: Logically Collective on Vec
2076: Input Parameter:
2077: . x - vector
2079: Output Parameters:
2080: + xx_v - the Fortran90 pointer to the array
2081: - ierr - error code
2083: Example of Usage:
2084: .vb
2085: #include <petsc/finclude/petscvec.h>
2086: use petscvec
2088: PetscScalar, pointer :: xx_v(:)
2089: ....
2090: call VecGetArrayReadF90(x,xx_v,ierr)
2091: a = xx_v(3)
2092: call VecRestoreArrayReadF90(x,xx_v,ierr)
2093: .ve
2095: If you intend to write entries into the array you must use VecGetArrayF90().
2097: Level: beginner
2099: .seealso: VecRestoreArrayReadF90(), VecGetArray(), VecRestoreArray(), VecGetArrayRead(), VecRestoreArrayRead(), VecGetArrayF90(), UsingFortran
2101: M*/
2103: /*MC
2104: VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2105: VecGetArrayReadF90().
2107: Synopsis:
2108: VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2110: Logically Collective on Vec
2112: Input Parameters:
2113: + x - vector
2114: - xx_v - the Fortran90 pointer to the array
2116: Output Parameter:
2117: . ierr - error code
2119: Example of Usage:
2120: .vb
2121: #include <petsc/finclude/petscvec.h>
2122: use petscvec
2124: PetscScalar, pointer :: xx_v(:)
2125: ....
2126: call VecGetArrayReadF90(x,xx_v,ierr)
2127: a = xx_v(3)
2128: call VecRestoreArrayReadF90(x,xx_v,ierr)
2129: .ve
2131: Level: beginner
2133: .seealso: VecGetArrayReadF90(), VecGetArray(), VecRestoreArray(), VecGetArrayRead(), VecRestoreArrayRead(),UsingFortran, VecRestoreArrayF90()
2135: M*/
2137: /*@C
2138: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2139: processor's portion of the vector data. You MUST call VecRestoreArray2d()
2140: when you no longer need access to the array.
2142: Logically Collective
2144: Input Parameter:
2145: + x - the vector
2146: . m - first dimension of two dimensional array
2147: . n - second dimension of two dimensional array
2148: . mstart - first index you will use in first coordinate direction (often 0)
2149: - nstart - first index in the second coordinate direction (often 0)
2151: Output Parameter:
2152: . a - location to put pointer to the array
2154: Level: developer
2156: Notes:
2157: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2158: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2159: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2160: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
2162: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2164: Concepts: vector^accessing local values as 2d array
2166: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2167: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2168: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2169: @*/
2170: PetscErrorCode VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2171: {
2173: PetscInt i,N;
2174: PetscScalar *aa;
2180: VecGetLocalSize(x,&N);
2181: if (m*n != N) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 2d array dimensions %D by %D",N,m,n);
2182: VecGetArray(x,&aa);
2184: PetscMalloc1(m,a);
2185: for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
2186: *a -= mstart;
2187: return(0);
2188: }
2190: /*@C
2191: VecRestoreArray2d - Restores a vector after VecGetArray2d() has been called.
2193: Logically Collective
2195: Input Parameters:
2196: + x - the vector
2197: . m - first dimension of two dimensional array
2198: . n - second dimension of the two dimensional array
2199: . mstart - first index you will use in first coordinate direction (often 0)
2200: . nstart - first index in the second coordinate direction (often 0)
2201: - a - location of pointer to array obtained from VecGetArray2d()
2203: Level: developer
2205: Notes:
2206: For regular PETSc vectors this routine does not involve any copies. For
2207: any special vectors that do not store local vector data in a contiguous
2208: array, this routine will copy the data back into the underlying
2209: vector data structure from the array obtained with VecGetArray().
2211: This routine actually zeros out the a pointer.
2213: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2214: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2215: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2216: @*/
2217: PetscErrorCode VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2218: {
2220: void *dummy;
2226: dummy = (void*)(*a + mstart);
2227: PetscFree(dummy);
2228: VecRestoreArray(x,NULL);
2229: return(0);
2230: }
2232: /*@C
2233: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2234: processor's portion of the vector data. You MUST call VecRestoreArray1d()
2235: when you no longer need access to the array.
2237: Logically Collective
2239: Input Parameter:
2240: + x - the vector
2241: . m - first dimension of two dimensional array
2242: - mstart - first index you will use in first coordinate direction (often 0)
2244: Output Parameter:
2245: . a - location to put pointer to the array
2247: Level: developer
2249: Notes:
2250: For a vector obtained from DMCreateLocalVector() mstart are likely
2251: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2252: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
2254: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2256: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2257: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2258: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2259: @*/
2260: PetscErrorCode VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2261: {
2263: PetscInt N;
2269: VecGetLocalSize(x,&N);
2270: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2271: VecGetArray(x,a);
2272: *a -= mstart;
2273: return(0);
2274: }
2276: /*@C
2277: VecRestoreArray1d - Restores a vector after VecGetArray1d() has been called.
2279: Logically Collective
2281: Input Parameters:
2282: + x - the vector
2283: . m - first dimension of two dimensional array
2284: . mstart - first index you will use in first coordinate direction (often 0)
2285: - a - location of pointer to array obtained from VecGetArray21()
2287: Level: developer
2289: Notes:
2290: For regular PETSc vectors this routine does not involve any copies. For
2291: any special vectors that do not store local vector data in a contiguous
2292: array, this routine will copy the data back into the underlying
2293: vector data structure from the array obtained with VecGetArray1d().
2295: This routine actually zeros out the a pointer.
2297: Concepts: vector^accessing local values as 1d array
2299: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2300: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2301: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2302: @*/
2303: PetscErrorCode VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2304: {
2310: VecRestoreArray(x,NULL);
2311: return(0);
2312: }
2315: /*@C
2316: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
2317: processor's portion of the vector data. You MUST call VecRestoreArray3d()
2318: when you no longer need access to the array.
2320: Logically Collective
2322: Input Parameter:
2323: + x - the vector
2324: . m - first dimension of three dimensional array
2325: . n - second dimension of three dimensional array
2326: . p - third dimension of three dimensional array
2327: . mstart - first index you will use in first coordinate direction (often 0)
2328: . nstart - first index in the second coordinate direction (often 0)
2329: - pstart - first index in the third coordinate direction (often 0)
2331: Output Parameter:
2332: . a - location to put pointer to the array
2334: Level: developer
2336: Notes:
2337: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
2338: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2339: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2340: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
2342: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2344: Concepts: vector^accessing local values as 3d array
2346: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2347: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2348: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2349: @*/
2350: PetscErrorCode VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2351: {
2353: PetscInt i,N,j;
2354: PetscScalar *aa,**b;
2360: VecGetLocalSize(x,&N);
2361: if (m*n*p != N) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 3d array dimensions %D by %D by %D",N,m,n,p);
2362: VecGetArray(x,&aa);
2364: PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
2365: b = (PetscScalar**)((*a) + m);
2366: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2367: for (i=0; i<m; i++)
2368: for (j=0; j<n; j++)
2369: b[i*n+j] = aa + i*n*p + j*p - pstart;
2371: *a -= mstart;
2372: return(0);
2373: }
2375: /*@C
2376: VecRestoreArray3d - Restores a vector after VecGetArray3d() has been called.
2378: Logically Collective
2380: Input Parameters:
2381: + x - the vector
2382: . m - first dimension of three dimensional array
2383: . n - second dimension of the three dimensional array
2384: . p - third dimension of the three dimensional array
2385: . mstart - first index you will use in first coordinate direction (often 0)
2386: . nstart - first index in the second coordinate direction (often 0)
2387: . pstart - first index in the third coordinate direction (often 0)
2388: - a - location of pointer to array obtained from VecGetArray3d()
2390: Level: developer
2392: Notes:
2393: For regular PETSc vectors this routine does not involve any copies. For
2394: any special vectors that do not store local vector data in a contiguous
2395: array, this routine will copy the data back into the underlying
2396: vector data structure from the array obtained with VecGetArray().
2398: This routine actually zeros out the a pointer.
2400: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2401: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2402: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2403: @*/
2404: PetscErrorCode VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2405: {
2407: void *dummy;
2413: dummy = (void*)(*a + mstart);
2414: PetscFree(dummy);
2415: VecRestoreArray(x,NULL);
2416: return(0);
2417: }
2419: /*@C
2420: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
2421: processor's portion of the vector data. You MUST call VecRestoreArray4d()
2422: when you no longer need access to the array.
2424: Logically Collective
2426: Input Parameter:
2427: + x - the vector
2428: . m - first dimension of four dimensional array
2429: . n - second dimension of four dimensional array
2430: . p - third dimension of four dimensional array
2431: . q - fourth dimension of four dimensional array
2432: . mstart - first index you will use in first coordinate direction (often 0)
2433: . nstart - first index in the second coordinate direction (often 0)
2434: . pstart - first index in the third coordinate direction (often 0)
2435: - qstart - first index in the fourth coordinate direction (often 0)
2437: Output Parameter:
2438: . a - location to put pointer to the array
2440: Level: beginner
2442: Notes:
2443: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
2444: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2445: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2446: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
2448: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2450: Concepts: vector^accessing local values as 3d array
2452: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2453: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2454: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2455: @*/
2456: PetscErrorCode VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2457: {
2459: PetscInt i,N,j,k;
2460: PetscScalar *aa,***b,**c;
2466: VecGetLocalSize(x,&N);
2467: if (m*n*p*q != N) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 4d array dimensions %D by %D by %D by %D",N,m,n,p,q);
2468: VecGetArray(x,&aa);
2470: PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
2471: b = (PetscScalar***)((*a) + m);
2472: c = (PetscScalar**)(b + m*n);
2473: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2474: for (i=0; i<m; i++)
2475: for (j=0; j<n; j++)
2476: b[i*n+j] = c + i*n*p + j*p - pstart;
2477: for (i=0; i<m; i++)
2478: for (j=0; j<n; j++)
2479: for (k=0; k<p; k++)
2480: c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
2481: *a -= mstart;
2482: return(0);
2483: }
2485: /*@C
2486: VecRestoreArray4d - Restores a vector after VecGetArray3d() has been called.
2488: Logically Collective
2490: Input Parameters:
2491: + x - the vector
2492: . m - first dimension of four dimensional array
2493: . n - second dimension of the four dimensional array
2494: . p - third dimension of the four dimensional array
2495: . q - fourth dimension of the four dimensional array
2496: . mstart - first index you will use in first coordinate direction (often 0)
2497: . nstart - first index in the second coordinate direction (often 0)
2498: . pstart - first index in the third coordinate direction (often 0)
2499: . qstart - first index in the fourth coordinate direction (often 0)
2500: - a - location of pointer to array obtained from VecGetArray4d()
2502: Level: beginner
2504: Notes:
2505: For regular PETSc vectors this routine does not involve any copies. For
2506: any special vectors that do not store local vector data in a contiguous
2507: array, this routine will copy the data back into the underlying
2508: vector data structure from the array obtained with VecGetArray().
2510: This routine actually zeros out the a pointer.
2512: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2513: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2514: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2515: @*/
2516: PetscErrorCode VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2517: {
2519: void *dummy;
2525: dummy = (void*)(*a + mstart);
2526: PetscFree(dummy);
2527: VecRestoreArray(x,NULL);
2528: return(0);
2529: }
2531: /*@C
2532: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
2533: processor's portion of the vector data. You MUST call VecRestoreArray2dRead()
2534: when you no longer need access to the array.
2536: Logically Collective
2538: Input Parameter:
2539: + x - the vector
2540: . m - first dimension of two dimensional array
2541: . n - second dimension of two dimensional array
2542: . mstart - first index you will use in first coordinate direction (often 0)
2543: - nstart - first index in the second coordinate direction (often 0)
2545: Output Parameter:
2546: . a - location to put pointer to the array
2548: Level: developer
2550: Notes:
2551: For a vector obtained from DMCreateLocalVector() mstart and nstart are likely
2552: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2553: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2554: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
2556: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2558: Concepts: vector^accessing local values as 2d array
2560: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2561: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2562: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2563: @*/
2564: PetscErrorCode VecGetArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2565: {
2566: PetscErrorCode ierr;
2567: PetscInt i,N;
2568: const PetscScalar *aa;
2574: VecGetLocalSize(x,&N);
2575: if (m*n != N) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 2d array dimensions %D by %D",N,m,n);
2576: VecGetArrayRead(x,&aa);
2578: PetscMalloc1(m,a);
2579: for (i=0; i<m; i++) (*a)[i] = (PetscScalar*) aa + i*n - nstart;
2580: *a -= mstart;
2581: return(0);
2582: }
2584: /*@C
2585: VecRestoreArray2dRead - Restores a vector after VecGetArray2dRead() has been called.
2587: Logically Collective
2589: Input Parameters:
2590: + x - the vector
2591: . m - first dimension of two dimensional array
2592: . n - second dimension of the two dimensional array
2593: . mstart - first index you will use in first coordinate direction (often 0)
2594: . nstart - first index in the second coordinate direction (often 0)
2595: - a - location of pointer to array obtained from VecGetArray2d()
2597: Level: developer
2599: Notes:
2600: For regular PETSc vectors this routine does not involve any copies. For
2601: any special vectors that do not store local vector data in a contiguous
2602: array, this routine will copy the data back into the underlying
2603: vector data structure from the array obtained with VecGetArray().
2605: This routine actually zeros out the a pointer.
2607: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2608: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2609: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2610: @*/
2611: PetscErrorCode VecRestoreArray2dRead(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
2612: {
2614: void *dummy;
2620: dummy = (void*)(*a + mstart);
2621: PetscFree(dummy);
2622: VecRestoreArrayRead(x,NULL);
2623: return(0);
2624: }
2626: /*@C
2627: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
2628: processor's portion of the vector data. You MUST call VecRestoreArray1dRead()
2629: when you no longer need access to the array.
2631: Logically Collective
2633: Input Parameter:
2634: + x - the vector
2635: . m - first dimension of two dimensional array
2636: - mstart - first index you will use in first coordinate direction (often 0)
2638: Output Parameter:
2639: . a - location to put pointer to the array
2641: Level: developer
2643: Notes:
2644: For a vector obtained from DMCreateLocalVector() mstart are likely
2645: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2646: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners().
2648: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2650: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2651: VecRestoreArray2d(), DMDAVecGetArray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2652: VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2653: @*/
2654: PetscErrorCode VecGetArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2655: {
2657: PetscInt N;
2663: VecGetLocalSize(x,&N);
2664: if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
2665: VecGetArrayRead(x,(const PetscScalar**)a);
2666: *a -= mstart;
2667: return(0);
2668: }
2670: /*@C
2671: VecRestoreArray1dRead - Restores a vector after VecGetArray1dRead() has been called.
2673: Logically Collective
2675: Input Parameters:
2676: + x - the vector
2677: . m - first dimension of two dimensional array
2678: . mstart - first index you will use in first coordinate direction (often 0)
2679: - a - location of pointer to array obtained from VecGetArray21()
2681: Level: developer
2683: Notes:
2684: For regular PETSc vectors this routine does not involve any copies. For
2685: any special vectors that do not store local vector data in a contiguous
2686: array, this routine will copy the data back into the underlying
2687: vector data structure from the array obtained with VecGetArray1dRead().
2689: This routine actually zeros out the a pointer.
2691: Concepts: vector^accessing local values as 1d array
2693: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2694: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2695: VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
2696: @*/
2697: PetscErrorCode VecRestoreArray1dRead(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
2698: {
2704: VecRestoreArrayRead(x,NULL);
2705: return(0);
2706: }
2709: /*@C
2710: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
2711: processor's portion of the vector data. You MUST call VecRestoreArray3dRead()
2712: when you no longer need access to the array.
2714: Logically Collective
2716: Input Parameter:
2717: + x - the vector
2718: . m - first dimension of three dimensional array
2719: . n - second dimension of three dimensional array
2720: . p - third dimension of three dimensional array
2721: . mstart - first index you will use in first coordinate direction (often 0)
2722: . nstart - first index in the second coordinate direction (often 0)
2723: - pstart - first index in the third coordinate direction (often 0)
2725: Output Parameter:
2726: . a - location to put pointer to the array
2728: Level: developer
2730: Notes:
2731: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
2732: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2733: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2734: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3dRead().
2736: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2738: Concepts: vector^accessing local values as 3d array
2740: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2741: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2742: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2743: @*/
2744: PetscErrorCode VecGetArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2745: {
2746: PetscErrorCode ierr;
2747: PetscInt i,N,j;
2748: const PetscScalar *aa;
2749: PetscScalar **b;
2755: VecGetLocalSize(x,&N);
2756: if (m*n*p != N) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 3d array dimensions %D by %D by %D",N,m,n,p);
2757: VecGetArrayRead(x,&aa);
2759: PetscMalloc1(m*sizeof(PetscScalar**)+m*n,a);
2760: b = (PetscScalar**)((*a) + m);
2761: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2762: for (i=0; i<m; i++)
2763: for (j=0; j<n; j++)
2764: b[i*n+j] = (PetscScalar *)aa + i*n*p + j*p - pstart;
2766: *a -= mstart;
2767: return(0);
2768: }
2770: /*@C
2771: VecRestoreArray3dRead - Restores a vector after VecGetArray3dRead() has been called.
2773: Logically Collective
2775: Input Parameters:
2776: + x - the vector
2777: . m - first dimension of three dimensional array
2778: . n - second dimension of the three dimensional array
2779: . p - third dimension of the three dimensional array
2780: . mstart - first index you will use in first coordinate direction (often 0)
2781: . nstart - first index in the second coordinate direction (often 0)
2782: . pstart - first index in the third coordinate direction (often 0)
2783: - a - location of pointer to array obtained from VecGetArray3dRead()
2785: Level: developer
2787: Notes:
2788: For regular PETSc vectors this routine does not involve any copies. For
2789: any special vectors that do not store local vector data in a contiguous
2790: array, this routine will copy the data back into the underlying
2791: vector data structure from the array obtained with VecGetArray().
2793: This routine actually zeros out the a pointer.
2795: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2796: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2797: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2798: @*/
2799: PetscErrorCode VecRestoreArray3dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
2800: {
2802: void *dummy;
2808: dummy = (void*)(*a + mstart);
2809: PetscFree(dummy);
2810: VecRestoreArrayRead(x,NULL);
2811: return(0);
2812: }
2814: /*@C
2815: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
2816: processor's portion of the vector data. You MUST call VecRestoreArray4dRead()
2817: when you no longer need access to the array.
2819: Logically Collective
2821: Input Parameter:
2822: + x - the vector
2823: . m - first dimension of four dimensional array
2824: . n - second dimension of four dimensional array
2825: . p - third dimension of four dimensional array
2826: . q - fourth dimension of four dimensional array
2827: . mstart - first index you will use in first coordinate direction (often 0)
2828: . nstart - first index in the second coordinate direction (often 0)
2829: . pstart - first index in the third coordinate direction (often 0)
2830: - qstart - first index in the fourth coordinate direction (often 0)
2832: Output Parameter:
2833: . a - location to put pointer to the array
2835: Level: beginner
2837: Notes:
2838: For a vector obtained from DMCreateLocalVector() mstart, nstart, and pstart are likely
2839: obtained from the corner indices obtained from DMDAGetGhostCorners() while for
2840: DMCreateGlobalVector() they are the corner indices from DMDAGetCorners(). In both cases
2841: the arguments from DMDAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
2843: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2845: Concepts: vector^accessing local values as 3d array
2847: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
2848: VecRestoreArray2d(), DMDAVecGetarray(), DMDAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
2849: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
2850: @*/
2851: PetscErrorCode VecGetArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2852: {
2853: PetscErrorCode ierr;
2854: PetscInt i,N,j,k;
2855: const PetscScalar *aa;
2856: PetscScalar ***b,**c;
2862: VecGetLocalSize(x,&N);
2863: if (m*n*p*q != N) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 4d array dimensions %D by %D by %D by %D",N,m,n,p,q);
2864: VecGetArrayRead(x,&aa);
2866: PetscMalloc1(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p,a);
2867: b = (PetscScalar***)((*a) + m);
2868: c = (PetscScalar**)(b + m*n);
2869: for (i=0; i<m; i++) (*a)[i] = b + i*n - nstart;
2870: for (i=0; i<m; i++)
2871: for (j=0; j<n; j++)
2872: b[i*n+j] = c + i*n*p + j*p - pstart;
2873: for (i=0; i<m; i++)
2874: for (j=0; j<n; j++)
2875: for (k=0; k<p; k++)
2876: c[i*n*p+j*p+k] = (PetscScalar*) aa + i*n*p*q + j*p*q + k*q - qstart;
2877: *a -= mstart;
2878: return(0);
2879: }
2881: /*@C
2882: VecRestoreArray4dRead - Restores a vector after VecGetArray3d() has been called.
2884: Logically Collective
2886: Input Parameters:
2887: + x - the vector
2888: . m - first dimension of four dimensional array
2889: . n - second dimension of the four dimensional array
2890: . p - third dimension of the four dimensional array
2891: . q - fourth dimension of the four dimensional array
2892: . mstart - first index you will use in first coordinate direction (often 0)
2893: . nstart - first index in the second coordinate direction (often 0)
2894: . pstart - first index in the third coordinate direction (often 0)
2895: . qstart - first index in the fourth coordinate direction (often 0)
2896: - a - location of pointer to array obtained from VecGetArray4dRead()
2898: Level: beginner
2900: Notes:
2901: For regular PETSc vectors this routine does not involve any copies. For
2902: any special vectors that do not store local vector data in a contiguous
2903: array, this routine will copy the data back into the underlying
2904: vector data structure from the array obtained with VecGetArray().
2906: This routine actually zeros out the a pointer.
2908: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
2909: VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DMDAVecGetArray(), DMDAVecRestoreArray()
2910: VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
2911: @*/
2912: PetscErrorCode VecRestoreArray4dRead(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
2913: {
2915: void *dummy;
2921: dummy = (void*)(*a + mstart);
2922: PetscFree(dummy);
2923: VecRestoreArrayRead(x,NULL);
2924: return(0);
2925: }
2927: #if defined(PETSC_USE_DEBUG)
2929: /*@
2930: VecLockGet - Gets the current lock status of a vector
2932: Logically Collective on Vec
2934: Input Parameter:
2935: . x - the vector
2937: Output Parameter:
2938: . state - greater than zero indicates the vector is locked for read; less then zero indicates the vector is
2939: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
2941: Level: beginner
2943: Concepts: vector^accessing local values
2945: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop()
2946: @*/
2947: PetscErrorCode VecLockGet(Vec x,PetscInt *state)
2948: {
2951: *state = x->lock;
2952: return(0);
2953: }
2955: /*@
2956: VecLockReadPush - Pushes a read-only lock on a vector to prevent it from writing
2958: Logically Collective on Vec
2960: Input Parameter:
2961: . x - the vector
2963: Notes:
2964: If this is set then calls to VecGetArray() or VecSetValues() or any other routines that change the vectors values will fail.
2966: The call can be nested, i.e., called multiple times on the same vector, but each VecLockReadPush(x) has to have one matching
2967: VecLockReadPop(x), which removes the latest read-only lock.
2969: Level: beginner
2971: Concepts: vector^accessing local values
2973: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPop(), VecLockGet()
2974: @*/
2975: PetscErrorCode VecLockReadPush(Vec x)
2976: {
2979: if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector is already locked for exclusive write access but you want to read it");
2980: x->lock++;
2981: return(0);
2982: }
2984: /*@
2985: VecLockReadPop - Pops a read-only lock from a vector
2987: Logically Collective on Vec
2989: Input Parameter:
2990: . x - the vector
2992: Level: beginner
2994: Concepts: vector^accessing local values
2996: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockGet()
2997: @*/
2998: PetscErrorCode VecLockReadPop(Vec x)
2999: {
3002: x->lock--;
3003: if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector has been unlocked from read-only access too many times");
3004: return(0);
3005: }
3007: /*@C
3008: VecLockWriteSet_Private - Lock or unlock a vector for exclusive read/write access
3010: Logically Collective on Vec
3012: Input Parameter:
3013: + x - the vector
3014: - flg - PETSC_TRUE to lock the vector for writing; PETSC_FALSE to unlock it.
3016: Notes:
3017: The function is usefull in split-phase computations, which usually have a begin phase and an end phase.
3018: One can call VecLockWriteSet_Private(x,PETSC_TRUE) in the begin phase to lock a vector for exclusive
3019: access, and call VecLockWriteSet_Private(x,PETSC_FALSE) in the end phase to unlock the vector from exclusive
3020: access. In this way, one is ensured no other operations can access the vector in between. The code may like
3023: VecGetArray(x,&xdata); // begin phase
3024: VecLockWriteSet_Private(v,PETSC_TRUE);
3026: Other operations, which can not acceess x anymore (they can access xdata, of course)
3028: VecRestoreArray(x,&vdata); // end phase
3029: VecLockWriteSet_Private(v,PETSC_FALSE);
3031: The call can not be nested on the same vector, in other words, one can not call VecLockWriteSet_Private(x,PETSC_TRUE)
3032: again before calling VecLockWriteSet_Private(v,PETSC_FALSE).
3034: Level: beginner
3036: Concepts: vector^accessing local values
3038: .seealso: VecRestoreArray(), VecGetArrayRead(), VecLockReadPush(), VecLockReadPop(), VecLockGet()
3039: @*/
3040: PetscErrorCode VecLockWriteSet_Private(Vec x,PetscBool flg)
3041: {
3044: if (flg) {
3045: if (x->lock > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector is already locked for read-only access but you want to write it");
3046: else if (x->lock < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector is already locked for exclusive write access but you want to write it");
3047: else x->lock = -1;
3048: } else {
3049: if (x->lock != -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Vector is not locked for exclusive write access but you want to unlock it from that");
3050: x->lock = 0;
3051: }
3052: return(0);
3053: }
3055: /*@
3056: VecLockPush - Pushes a read-only lock on a vector to prevent it from writing
3058: Level: deprecated
3060: Concepts: vector^accessing local values
3062: .seealso: VecLockReadPush()
3063: @*/
3064: PetscErrorCode VecLockPush(Vec x)
3065: {
3068: VecLockReadPush(x);
3069: return(0);
3070: }
3072: /*@
3073: VecLockPop - Pops a read-only lock from a vector
3075: Level: deprecated
3077: Concepts: vector^accessing local values
3079: .seealso: VecLockReadPop()
3080: @*/
3081: PetscErrorCode VecLockPop(Vec x)
3082: {
3085: VecLockReadPop(x);
3086: return(0);
3087: }
3089: #endif