Actual source code: rvector.c
1: /*
2: Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
3: These are the vector functions the user calls.
4: */
5: #include <petsc/private/vecimpl.h>
7: PetscInt VecGetSubVectorSavedStateId = -1;
9: #if PetscDefined(USE_DEBUG)
10: // this is a no-op '0' macro in optimized builds
11: PetscErrorCode VecValidValues_Internal(Vec vec, PetscInt argnum, PetscBool begin)
12: {
13: PetscFunctionBegin;
14: if (vec->petscnative || vec->ops->getarray) {
15: PetscInt n;
16: const PetscScalar *x;
17: PetscOffloadMask mask;
19: PetscCall(VecGetOffloadMask(vec, &mask));
20: if (!PetscOffloadHost(mask)) PetscFunctionReturn(PETSC_SUCCESS);
21: PetscCall(VecGetLocalSize(vec, &n));
22: PetscCall(VecGetArrayRead(vec, &x));
23: for (PetscInt i = 0; i < n; i++) {
24: if (begin) {
25: PetscCheck(!PetscIsInfOrNanScalar(x[i]), PETSC_COMM_SELF, PETSC_ERR_FP, "Vec entry at local location %" PetscInt_FMT " is not-a-number or infinite at beginning of function: Parameter number %" PetscInt_FMT, i, argnum);
26: } else {
27: PetscCheck(!PetscIsInfOrNanScalar(x[i]), PETSC_COMM_SELF, PETSC_ERR_FP, "Vec entry at local location %" PetscInt_FMT " is not-a-number or infinite at end of function: Parameter number %" PetscInt_FMT, i, argnum);
28: }
29: }
30: PetscCall(VecRestoreArrayRead(vec, &x));
31: }
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
34: #endif
36: /*@
37: VecMaxPointwiseDivide - Computes the maximum of the componentwise division `max = max_i abs(x[i]/y[i])`.
39: Logically Collective
41: Input Parameters:
42: + x - the numerators
43: - y - the denominators
45: Output Parameter:
46: . max - the result
48: Level: advanced
50: Notes:
51: `x` and `y` may be the same vector
53: if a particular `y[i]` is zero, it is treated as 1 in the above formula
55: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`
56: @*/
57: PetscErrorCode VecMaxPointwiseDivide(Vec x, Vec y, PetscReal *max)
58: {
59: PetscFunctionBegin;
62: PetscAssertPointer(max, 3);
65: PetscCheckSameTypeAndComm(x, 1, y, 2);
66: VecCheckSameSize(x, 1, y, 2);
67: VecCheckAssembled(x);
68: VecCheckAssembled(y);
69: PetscCall(VecLockReadPush(x));
70: PetscCall(VecLockReadPush(y));
71: PetscUseTypeMethod(x, maxpointwisedivide, y, max);
72: PetscCall(VecLockReadPop(x));
73: PetscCall(VecLockReadPop(y));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: /*@
78: VecDot - Computes the vector dot product.
80: Collective
82: Input Parameters:
83: + x - first vector
84: - y - second vector
86: Output Parameter:
87: . val - the dot product
89: Level: intermediate
91: Notes for Users of Complex Numbers:
92: For complex vectors, `VecDot()` computes
93: $ val = (x,y) = y^H x,
94: where y^H denotes the conjugate transpose of y. Note that this corresponds to the usual "mathematicians" complex
95: inner product where the SECOND argument gets the complex conjugate. Since the `BLASdot()` complex conjugates the first
96: first argument we call the `BLASdot()` with the arguments reversed.
98: Use `VecTDot()` for the indefinite form
99: $ val = (x,y) = y^T x,
100: where y^T denotes the transpose of y.
102: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDotRealPart()`
103: @*/
104: PetscErrorCode VecDot(Vec x, Vec y, PetscScalar *val)
105: {
106: PetscFunctionBegin;
109: PetscAssertPointer(val, 3);
112: PetscCheckSameTypeAndComm(x, 1, y, 2);
113: VecCheckSameSize(x, 1, y, 2);
114: VecCheckAssembled(x);
115: VecCheckAssembled(y);
117: PetscCall(VecLockReadPush(x));
118: PetscCall(VecLockReadPush(y));
119: PetscCall(PetscLogEventBegin(VEC_Dot, x, y, 0, 0));
120: PetscUseTypeMethod(x, dot, y, val);
121: PetscCall(PetscLogEventEnd(VEC_Dot, x, y, 0, 0));
122: PetscCall(VecLockReadPop(x));
123: PetscCall(VecLockReadPop(y));
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: /*@
128: VecDotRealPart - Computes the real part of the vector dot product.
130: Collective
132: Input Parameters:
133: + x - first vector
134: - y - second vector
136: Output Parameter:
137: . val - the real part of the dot product;
139: Level: intermediate
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 Notes:
150: This is not currently optimized to compute only the real part of the dot product.
152: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDot()`, `VecDotNorm2()`
153: @*/
154: PetscErrorCode VecDotRealPart(Vec x, Vec y, PetscReal *val)
155: {
156: PetscScalar fdot;
158: PetscFunctionBegin;
159: PetscCall(VecDot(x, y, &fdot));
160: *val = PetscRealPart(fdot);
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: /*@
165: VecNorm - Computes the vector norm.
167: Collective
169: Input Parameters:
170: + x - the vector
171: - type - the type of the norm requested
173: Output Parameter:
174: . val - the norm
176: Level: intermediate
178: Notes:
179: See `NormType` for descriptions of each norm.
181: For complex numbers `NORM_1` will return the traditional 1 norm of the 2 norm of the complex
182: numbers; that is the 1 norm of the absolute values of the complex entries. In PETSc 3.6 and
183: earlier releases it returned the 1 norm of the 1 norm of the complex entries (what is
184: returned by the BLAS routine `asum()`). Both are valid norms but most people expect the former.
186: This routine stashes the computed norm value, repeated calls before the vector entries are
187: changed are then rapid since the precomputed value is immediately available. Certain vector
188: operations such as `VecSet()` store the norms so the value is immediately available and does
189: not need to be explicitly computed. `VecScale()` updates any stashed norm values, thus calls
190: after `VecScale()` do not need to explicitly recompute the norm.
192: .seealso: [](ch_vectors), `Vec`, `NormType`, `VecDot()`, `VecTDot()`, `VecDotBegin()`, `VecDotEnd()`, `VecNormAvailable()`,
193: `VecNormBegin()`, `VecNormEnd()`, `NormType()`
194: @*/
195: PetscErrorCode VecNorm(Vec x, NormType type, PetscReal *val)
196: {
197: PetscBool flg = PETSC_TRUE;
199: PetscFunctionBegin;
202: VecCheckAssembled(x);
204: PetscAssertPointer(val, 3);
206: PetscCall(VecNormAvailable(x, type, &flg, val));
207: // check that all MPI processes call this routine together and have same availability
208: if (PetscDefined(USE_DEBUG)) {
209: PetscMPIInt b0 = (PetscMPIInt)flg, b1[2], b2[2];
210: b1[0] = -b0;
211: b1[1] = b0;
212: PetscCallMPI(MPIU_Allreduce(b1, b2, 2, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
213: PetscCheck(-b2[0] == b2[1], PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_WRONGSTATE, "Some MPI processes have cached %s norm, others do not. This may happen when some MPI processes call VecGetArray() and some others do not.", NormTypes[type]);
214: if (flg) {
215: PetscReal b1[2], b2[2];
216: b1[0] = -(*val);
217: b1[1] = *val;
218: PetscCallMPI(MPIU_Allreduce(b1, b2, 2, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)x)));
219: PetscCheck((PetscIsNanReal(b2[0]) && PetscIsNanReal(b2[1])) || (-b2[0] == b2[1]), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Difference in cached %s norms: local %g", NormTypes[type], (double)*val);
220: PetscCheck(-b2[0] == b2[1], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Difference in cached %s norms: local %g", NormTypes[type], (double)*val);
221: }
222: }
223: if (flg) PetscFunctionReturn(PETSC_SUCCESS);
225: PetscCall(VecLockReadPush(x));
226: PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
227: PetscUseTypeMethod(x, norm, type, val);
228: PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
229: PetscCall(VecLockReadPop(x));
231: if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: /*@
236: VecNormAvailable - Returns the vector norm if it is already known. That is, it has been previously computed and cached in the vector
238: Not Collective
240: Input Parameters:
241: + x - the vector
242: - type - one of `NORM_1` (sum_i |x[i]|), `NORM_2` sqrt(sum_i (x[i])^2), `NORM_INFINITY` max_i |x[i]|. Also available
243: `NORM_1_AND_2`, which computes both norms and stores them
244: in a two element array.
246: Output Parameters:
247: + available - `PETSC_TRUE` if the val returned is valid
248: - val - the norm
250: Level: intermediate
252: Developer Notes:
253: `PETSC_HAVE_SLOW_BLAS_NORM2` will cause a C (loop unrolled) version of the norm to be used, rather
254: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
255: (as opposed to vendor provided) because the FORTRAN BLAS `NRM2()` routine is very slow.
257: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
258: `VecNormBegin()`, `VecNormEnd()`
259: @*/
260: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
261: {
262: PetscFunctionBegin;
265: PetscAssertPointer(available, 3);
266: PetscAssertPointer(val, 4);
268: if (type == NORM_1_AND_2) {
269: *available = PETSC_FALSE;
270: } else {
271: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
272: }
273: PetscFunctionReturn(PETSC_SUCCESS);
274: }
276: /*@
277: VecNormalize - Normalizes a vector by its 2-norm.
279: Collective
281: Input Parameter:
282: . x - the vector
284: Output Parameter:
285: . val - the vector norm before normalization. May be `NULL` if the value is not needed.
287: Level: intermediate
289: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
290: @*/
291: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
292: {
293: PetscReal norm;
295: PetscFunctionBegin;
298: PetscCall(VecSetErrorIfLocked(x, 1));
299: if (val) PetscAssertPointer(val, 2);
300: PetscCall(PetscLogEventBegin(VEC_Normalize, x, 0, 0, 0));
301: PetscCall(VecNorm(x, NORM_2, &norm));
302: if (norm == 0.0) PetscCall(PetscInfo(x, "Vector of zero norm can not be normalized; Returning only the zero norm\n"));
303: else if (PetscIsInfOrNanReal(norm)) PetscCall(PetscInfo(x, "Vector with Inf or Nan norm can not be normalized; Returning only the norm\n"));
304: else {
305: PetscScalar s = 1.0 / norm;
306: PetscCall(VecScale(x, s));
307: }
308: PetscCall(PetscLogEventEnd(VEC_Normalize, x, 0, 0, 0));
309: if (val) *val = norm;
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: /*@
314: VecMax - Determines the vector component with maximum real part and its location.
316: Collective
318: Input Parameter:
319: . x - the vector
321: Output Parameters:
322: + p - the index of `val` (pass `NULL` if you don't want this) in the vector
323: - val - the maximum component
325: Level: intermediate
327: Notes:
328: Returns the value `PETSC_MIN_REAL` and negative `p` if the vector is of length 0.
330: Returns the smallest index with the maximum value
332: Developer Note:
333: The Nag Fortran compiler does not like the symbol name VecMax
335: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
336: @*/
337: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
338: {
339: PetscFunctionBegin;
342: VecCheckAssembled(x);
343: if (p) PetscAssertPointer(p, 2);
344: PetscAssertPointer(val, 3);
345: PetscCall(VecLockReadPush(x));
346: PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
347: PetscUseTypeMethod(x, max, p, val);
348: PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
349: PetscCall(VecLockReadPop(x));
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }
353: /*@
354: VecMin - Determines the vector component with minimum real part and its location.
356: Collective
358: Input Parameter:
359: . x - the vector
361: Output Parameters:
362: + p - the index of `val` (pass `NULL` if you don't want this location) in the vector
363: - val - the minimum component
365: Level: intermediate
367: Notes:
368: Returns the value `PETSC_MAX_REAL` and negative `p` if the vector is of length 0.
370: This returns the smallest index with the minimum value
372: Developer Note:
373: The Nag Fortran compiler does not like the symbol name VecMin
375: .seealso: [](ch_vectors), `Vec`, `VecMax()`
376: @*/
377: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
378: {
379: PetscFunctionBegin;
382: VecCheckAssembled(x);
383: if (p) PetscAssertPointer(p, 2);
384: PetscAssertPointer(val, 3);
385: PetscCall(VecLockReadPush(x));
386: PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
387: PetscUseTypeMethod(x, min, p, val);
388: PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
389: PetscCall(VecLockReadPop(x));
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: /*@
394: VecTDot - Computes an indefinite vector dot product. That is, this
395: routine does NOT use the complex conjugate.
397: Collective
399: Input Parameters:
400: + x - first vector
401: - y - second vector
403: Output Parameter:
404: . val - the dot product
406: Level: intermediate
408: Notes for Users of Complex Numbers:
409: For complex vectors, `VecTDot()` computes the indefinite form
410: $ val = (x,y) = y^T x,
411: where y^T denotes the transpose of y.
413: Use `VecDot()` for the inner product
414: $ val = (x,y) = y^H x,
415: where y^H denotes the conjugate transpose of y.
417: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
418: @*/
419: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
420: {
421: PetscFunctionBegin;
424: PetscAssertPointer(val, 3);
427: PetscCheckSameTypeAndComm(x, 1, y, 2);
428: VecCheckSameSize(x, 1, y, 2);
429: VecCheckAssembled(x);
430: VecCheckAssembled(y);
432: PetscCall(VecLockReadPush(x));
433: PetscCall(VecLockReadPush(y));
434: PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
435: PetscUseTypeMethod(x, tdot, y, val);
436: PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
437: PetscCall(VecLockReadPop(x));
438: PetscCall(VecLockReadPop(y));
439: PetscFunctionReturn(PETSC_SUCCESS);
440: }
442: PetscErrorCode VecScaleAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
443: {
444: PetscReal norms[4];
445: PetscBool flgs[4];
446: PetscScalar one = 1.0;
448: PetscFunctionBegin;
451: VecCheckAssembled(x);
452: PetscCall(VecSetErrorIfLocked(x, 1));
454: if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);
456: /* get current stashed norms */
457: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
459: PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
460: VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
461: PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));
463: PetscCall(PetscObjectStateIncrease((PetscObject)x));
464: /* put the scaled stashed norms back into the Vec */
465: for (PetscInt i = 0; i < 4; i++) {
466: PetscReal ar = PetscAbsScalar(alpha);
467: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
468: }
469: PetscFunctionReturn(PETSC_SUCCESS);
470: }
472: /*@
473: VecScale - Scales a vector.
475: Logically Collective
477: Input Parameters:
478: + x - the vector
479: - alpha - the scalar
481: Level: intermediate
483: Note:
484: For a vector with n components, `VecScale()` computes x[i] = alpha * x[i], for i=1,...,n.
486: .seealso: [](ch_vectors), `Vec`, `VecSet()`
487: @*/
488: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
489: {
490: PetscFunctionBegin;
491: PetscCall(VecScaleAsync_Private(x, alpha, NULL));
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
496: {
497: PetscFunctionBegin;
500: VecCheckAssembled(x);
502: PetscCall(VecSetErrorIfLocked(x, 1));
504: if (alpha == 0) {
505: PetscReal norm;
506: PetscBool set;
508: PetscCall(VecNormAvailable(x, NORM_2, &set, &norm));
509: if (set == PETSC_TRUE && norm == 0) PetscFunctionReturn(PETSC_SUCCESS);
510: }
511: PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
512: VecMethodDispatch(x, dctx, VecAsyncFnName(Set), set, (Vec, PetscScalar, PetscDeviceContext), alpha);
513: PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
514: PetscCall(PetscObjectStateIncrease((PetscObject)x));
516: /* norms can be simply set (if |alpha|*N not too large) */
517: {
518: PetscReal val = PetscAbsScalar(alpha);
519: const PetscInt N = x->map->N;
521: if (N == 0) {
522: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0));
523: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
524: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
525: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
526: } else if (val > PETSC_MAX_REAL / N) {
527: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
528: } else {
529: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
530: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
531: val *= PetscSqrtReal((PetscReal)N);
532: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
533: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
534: }
535: }
536: PetscFunctionReturn(PETSC_SUCCESS);
537: }
539: /*@
540: VecSet - Sets all components of a vector to a single scalar value.
542: Logically Collective
544: Input Parameters:
545: + x - the vector
546: - alpha - the scalar
548: Level: beginner
550: Notes:
551: For a vector of dimension n, `VecSet()` sets x[i] = alpha, for i=1,...,n,
552: so that all vector entries then equal the identical
553: scalar value, `alpha`. Use the more general routine
554: `VecSetValues()` to set different vector entries.
556: You CANNOT call this after you have called `VecSetValues()` but before you call
557: `VecAssemblyBegin()`
559: If `alpha` is zero and the norm of the vector is known to be zero then this skips the unneeded zeroing process
561: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
562: @*/
563: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
564: {
565: PetscFunctionBegin;
566: PetscCall(VecSetAsync_Private(x, alpha, NULL));
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: PetscErrorCode VecAXPYAsync_Private(Vec y, PetscScalar alpha, Vec x, PetscDeviceContext dctx)
571: {
572: PetscFunctionBegin;
577: PetscCheckSameTypeAndComm(x, 3, y, 1);
578: VecCheckSameSize(x, 3, y, 1);
579: VecCheckAssembled(x);
580: VecCheckAssembled(y);
582: if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
583: PetscCall(VecSetErrorIfLocked(y, 1));
584: if (x == y) {
585: PetscCall(VecScale(y, alpha + 1.0));
586: PetscFunctionReturn(PETSC_SUCCESS);
587: }
588: PetscCall(VecLockReadPush(x));
589: PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
590: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPY), axpy, (Vec, PetscScalar, Vec, PetscDeviceContext), alpha, x);
591: PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
592: PetscCall(VecLockReadPop(x));
593: PetscCall(PetscObjectStateIncrease((PetscObject)y));
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
596: /*@
597: VecAXPY - Computes `y = alpha x + y`.
599: Logically Collective
601: Input Parameters:
602: + alpha - the scalar
603: . x - vector scale by `alpha`
604: - y - vector accumulated into
606: Output Parameter:
607: . y - output vector
609: Level: intermediate
611: Notes:
612: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
613: .vb
614: VecAXPY(y,alpha,x) y = alpha x + y
615: VecAYPX(y,beta,x) y = x + beta y
616: VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
617: VecWAXPY(w,alpha,x,y) w = alpha x + y
618: VecAXPBYPCZ(z,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
619: VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
620: .ve
622: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
623: @*/
624: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
625: {
626: PetscFunctionBegin;
627: PetscCall(VecAXPYAsync_Private(y, alpha, x, NULL));
628: PetscFunctionReturn(PETSC_SUCCESS);
629: }
631: PetscErrorCode VecAYPXAsync_Private(Vec y, PetscScalar beta, Vec x, PetscDeviceContext dctx)
632: {
633: PetscFunctionBegin;
638: PetscCheckSameTypeAndComm(x, 3, y, 1);
639: VecCheckSameSize(x, 1, y, 3);
640: VecCheckAssembled(x);
641: VecCheckAssembled(y);
643: PetscCall(VecSetErrorIfLocked(y, 1));
644: if (x == y) {
645: PetscCall(VecScale(y, beta + 1.0));
646: PetscFunctionReturn(PETSC_SUCCESS);
647: }
648: PetscCall(VecLockReadPush(x));
649: if (beta == (PetscScalar)0.0) {
650: PetscCall(VecCopy(x, y));
651: } else {
652: PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
653: VecMethodDispatch(y, dctx, VecAsyncFnName(AYPX), aypx, (Vec, PetscScalar, Vec, PetscDeviceContext), beta, x);
654: PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
655: PetscCall(PetscObjectStateIncrease((PetscObject)y));
656: }
657: PetscCall(VecLockReadPop(x));
658: PetscFunctionReturn(PETSC_SUCCESS);
659: }
661: /*@
662: VecAYPX - Computes `y = x + beta y`.
664: Logically Collective
666: Input Parameters:
667: + beta - the scalar
668: . x - the unscaled vector
669: - y - the vector to be scaled
671: Output Parameter:
672: . y - output vector
674: Level: intermediate
676: Developer Notes:
677: The implementation is optimized for `beta` of -1.0, 0.0, and 1.0
679: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
680: @*/
681: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
682: {
683: PetscFunctionBegin;
684: PetscCall(VecAYPXAsync_Private(y, beta, x, NULL));
685: PetscFunctionReturn(PETSC_SUCCESS);
686: }
688: PetscErrorCode VecAXPBYAsync_Private(Vec y, PetscScalar alpha, PetscScalar beta, Vec x, PetscDeviceContext dctx)
689: {
690: PetscFunctionBegin;
695: PetscCheckSameTypeAndComm(x, 4, y, 1);
696: VecCheckSameSize(y, 1, x, 4);
697: VecCheckAssembled(x);
698: VecCheckAssembled(y);
701: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
702: if (x == y) {
703: PetscCall(VecScale(y, alpha + beta));
704: PetscFunctionReturn(PETSC_SUCCESS);
705: }
707: PetscCall(VecSetErrorIfLocked(y, 1));
708: PetscCall(VecLockReadPush(x));
709: PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
710: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPBY), axpby, (Vec, PetscScalar, PetscScalar, Vec, PetscDeviceContext), alpha, beta, x);
711: PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
712: PetscCall(PetscObjectStateIncrease((PetscObject)y));
713: PetscCall(VecLockReadPop(x));
714: PetscFunctionReturn(PETSC_SUCCESS);
715: }
717: /*@
718: VecAXPBY - Computes `y = alpha x + beta y`.
720: Logically Collective
722: Input Parameters:
723: + alpha - first scalar
724: . beta - second scalar
725: . x - the first scaled vector
726: - y - the second scaled vector
728: Output Parameter:
729: . y - output vector
731: Level: intermediate
733: Developer Notes:
734: The implementation is optimized for `alpha` and/or `beta` values of 0.0 and 1.0
736: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
737: @*/
738: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
739: {
740: PetscFunctionBegin;
741: PetscCall(VecAXPBYAsync_Private(y, alpha, beta, x, NULL));
742: PetscFunctionReturn(PETSC_SUCCESS);
743: }
745: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
746: {
747: PetscFunctionBegin;
754: PetscCheckSameTypeAndComm(x, 5, y, 6);
755: PetscCheckSameTypeAndComm(x, 5, z, 1);
756: VecCheckSameSize(x, 5, y, 6);
757: VecCheckSameSize(x, 5, z, 1);
758: PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
759: PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
760: VecCheckAssembled(x);
761: VecCheckAssembled(y);
762: VecCheckAssembled(z);
766: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
768: PetscCall(VecSetErrorIfLocked(z, 1));
769: PetscCall(VecLockReadPush(x));
770: PetscCall(VecLockReadPush(y));
771: PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
772: VecMethodDispatch(z, dctx, VecAsyncFnName(AXPBYPCZ), axpbypcz, (Vec, PetscScalar, PetscScalar, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, beta, gamma, x, y);
773: PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
774: PetscCall(PetscObjectStateIncrease((PetscObject)z));
775: PetscCall(VecLockReadPop(x));
776: PetscCall(VecLockReadPop(y));
777: PetscFunctionReturn(PETSC_SUCCESS);
778: }
779: /*@
780: VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`
782: Logically Collective
784: Input Parameters:
785: + alpha - first scalar
786: . beta - second scalar
787: . gamma - third scalar
788: . x - first vector
789: . y - second vector
790: - z - third vector
792: Output Parameter:
793: . z - output vector
795: Level: intermediate
797: Note:
798: `x`, `y` and `z` must be different vectors
800: Developer Notes:
801: The implementation is optimized for `alpha` of 1.0 and `gamma` of 1.0 or 0.0
803: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
804: @*/
805: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
806: {
807: PetscFunctionBegin;
808: PetscCall(VecAXPBYPCZAsync_Private(z, alpha, beta, gamma, x, y, NULL));
809: PetscFunctionReturn(PETSC_SUCCESS);
810: }
812: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
813: {
814: PetscFunctionBegin;
821: PetscCheckSameTypeAndComm(x, 3, y, 4);
822: PetscCheckSameTypeAndComm(y, 4, w, 1);
823: VecCheckSameSize(x, 3, y, 4);
824: VecCheckSameSize(x, 3, w, 1);
825: PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
826: PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
827: VecCheckAssembled(x);
828: VecCheckAssembled(y);
830: PetscCall(VecSetErrorIfLocked(w, 1));
832: PetscCall(VecLockReadPush(x));
833: PetscCall(VecLockReadPush(y));
834: if (alpha == (PetscScalar)0.0) {
835: PetscCall(VecCopyAsync_Private(y, w, dctx));
836: } else {
837: PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
838: VecMethodDispatch(w, dctx, VecAsyncFnName(WAXPY), waxpy, (Vec, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, x, y);
839: PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
840: PetscCall(PetscObjectStateIncrease((PetscObject)w));
841: }
842: PetscCall(VecLockReadPop(x));
843: PetscCall(VecLockReadPop(y));
844: PetscFunctionReturn(PETSC_SUCCESS);
845: }
847: /*@
848: VecWAXPY - Computes `w = alpha x + y`.
850: Logically Collective
852: Input Parameters:
853: + alpha - the scalar
854: . x - first vector, multiplied by `alpha`
855: - y - second vector
857: Output Parameter:
858: . w - the result
860: Level: intermediate
862: Note:
863: `w` cannot be either `x` or `y`, but `x` and `y` can be the same
865: Developer Notes:
866: The implementation is optimized for alpha of -1.0, 0.0, and 1.0
868: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
869: @*/
870: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
871: {
872: PetscFunctionBegin;
873: PetscCall(VecWAXPYAsync_Private(w, alpha, x, y, NULL));
874: PetscFunctionReturn(PETSC_SUCCESS);
875: }
877: /*@
878: VecSetValues - Inserts or adds values into certain locations of a vector.
880: Not Collective
882: Input Parameters:
883: + x - vector to insert in
884: . ni - number of elements to add
885: . ix - indices where to add
886: . y - array of values
887: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries
889: Level: beginner
891: Notes:
892: .vb
893: `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
894: .ve
896: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
897: options cannot be mixed without intervening calls to the assembly
898: routines.
900: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
901: MUST be called after all calls to `VecSetValues()` have been completed.
903: VecSetValues() uses 0-based indices in Fortran as well as in C.
905: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
906: negative indices may be passed in ix. These rows are
907: simply ignored. This allows easily inserting element load matrices
908: with homogeneous Dirichlet boundary conditions that you don't want represented
909: in the vector.
911: Fortran Note:
912: If any of `ix` and `y` are scalars pass them using, for example,
913: .vb
914: call VecSetValues(mat, one, [ix], [y], INSERT_VALUES, ierr)
915: .ve
917: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
918: `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
919: @*/
920: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
921: {
922: PetscFunctionBeginHot;
924: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
925: PetscAssertPointer(ix, 3);
926: PetscAssertPointer(y, 4);
929: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
930: PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
931: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
932: PetscCall(PetscObjectStateIncrease((PetscObject)x));
933: PetscFunctionReturn(PETSC_SUCCESS);
934: }
936: /*@
937: VecGetValues - Gets values from certain locations of a vector. Currently
938: can only get values on the same processor on which they are owned
940: Not Collective
942: Input Parameters:
943: + x - vector to get values from
944: . ni - number of elements to get
945: - ix - indices where to get them from (in global 1d numbering)
947: Output Parameter:
948: . y - array of values, must be passed in with a length of `ni`
950: Level: beginner
952: Notes:
953: The user provides the allocated array y; it is NOT allocated in this routine
955: `VecGetValues()` gets y[i] = x[ix[i]], for i=0,...,ni-1.
957: `VecAssemblyBegin()` and `VecAssemblyEnd()` MUST be called before calling this if `VecSetValues()` or related routine has been called
959: VecGetValues() uses 0-based indices in Fortran as well as in C.
961: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
962: negative indices may be passed in ix. These rows are
963: simply ignored.
965: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
966: @*/
967: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
968: {
969: PetscFunctionBegin;
971: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
972: PetscAssertPointer(ix, 3);
973: PetscAssertPointer(y, 4);
975: VecCheckAssembled(x);
976: PetscUseTypeMethod(x, getvalues, ni, ix, y);
977: PetscFunctionReturn(PETSC_SUCCESS);
978: }
980: /*@
981: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
983: Not Collective
985: Input Parameters:
986: + x - vector to insert in
987: . ni - number of blocks to add
988: . ix - indices where to add in block count, rather than element count
989: . y - array of values
990: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries
992: Level: intermediate
994: Notes:
995: `VecSetValuesBlocked()` sets x[bs*ix[i]+j] = y[bs*i+j],
996: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
998: Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
999: options cannot be mixed without intervening calls to the assembly
1000: routines.
1002: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1003: MUST be called after all calls to `VecSetValuesBlocked()` have been completed.
1005: `VecSetValuesBlocked()` uses 0-based indices in Fortran as well as in C.
1007: Negative indices may be passed in ix, these rows are
1008: simply ignored. This allows easily inserting element load matrices
1009: with homogeneous Dirichlet boundary conditions that you don't want represented
1010: in the vector.
1012: Fortran Note:
1013: If any of `ix` and `y` are scalars pass them using, for example,
1014: .vb
1015: call VecSetValuesBlocked(mat, one, [ix], [y], INSERT_VALUES, ierr)
1016: .ve
1018: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
1019: `VecSetValues()`
1020: @*/
1021: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1022: {
1023: PetscFunctionBeginHot;
1025: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1026: PetscAssertPointer(ix, 3);
1027: PetscAssertPointer(y, 4);
1030: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1031: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1032: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1033: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1034: PetscFunctionReturn(PETSC_SUCCESS);
1035: }
1037: /*@
1038: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1039: using a local ordering of the nodes.
1041: Not Collective
1043: Input Parameters:
1044: + x - vector to insert in
1045: . ni - number of elements to add
1046: . ix - indices where to add
1047: . y - array of values
1048: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1050: Level: intermediate
1052: Notes:
1053: `VecSetValuesLocal()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
1055: Calls to `VecSetValuesLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1056: options cannot be mixed without intervening calls to the assembly
1057: routines.
1059: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1060: MUST be called after all calls to `VecSetValuesLocal()` have been completed.
1062: `VecSetValuesLocal()` uses 0-based indices in Fortran as well as in C.
1064: Fortran Note:
1065: If any of `ix` and `y` are scalars pass them using, for example,
1066: .vb
1067: call VecSetValuesLocal(mat, one, [ix], [y], INSERT_VALUES, ierr)
1068: .ve
1070: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1071: `VecSetValuesBlockedLocal()`
1072: @*/
1073: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1074: {
1075: PetscInt lixp[128], *lix = lixp;
1077: PetscFunctionBeginHot;
1079: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1080: PetscAssertPointer(ix, 3);
1081: PetscAssertPointer(y, 4);
1084: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1085: if (!x->ops->setvalueslocal) {
1086: if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
1087: if (x->map->mapping) {
1088: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1089: PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1090: PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1091: if (ni > 128) PetscCall(PetscFree(lix));
1092: } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1093: } else PetscUseTypeMethod(x, setvalueslocal, ni, ix, y, iora);
1094: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1095: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1096: PetscFunctionReturn(PETSC_SUCCESS);
1097: }
1099: /*@
1100: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1101: using a local ordering of the nodes.
1103: Not Collective
1105: Input Parameters:
1106: + x - vector to insert in
1107: . ni - number of blocks to add
1108: . ix - indices where to add in block count, not element count
1109: . y - array of values
1110: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1112: Level: intermediate
1114: Notes:
1115: `VecSetValuesBlockedLocal()` sets x[bs*ix[i]+j] = y[bs*i+j],
1116: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with `VecSetBlockSize()`.
1118: Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1119: options cannot be mixed without intervening calls to the assembly
1120: routines.
1122: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1123: MUST be called after all calls to `VecSetValuesBlockedLocal()` have been completed.
1125: `VecSetValuesBlockedLocal()` uses 0-based indices in Fortran as well as in C.
1127: Fortran Note:
1128: If any of `ix` and `y` are scalars pass them using, for example,
1129: .vb
1130: call VecSetValuesBlockedLocal(mat, one, [ix], [y], INSERT_VALUES, ierr)
1131: .ve
1133: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1134: `VecSetLocalToGlobalMapping()`
1135: @*/
1136: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1137: {
1138: PetscInt lixp[128], *lix = lixp;
1140: PetscFunctionBeginHot;
1142: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1143: PetscAssertPointer(ix, 3);
1144: PetscAssertPointer(y, 4);
1146: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1147: if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
1148: if (x->map->mapping) {
1149: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1150: PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1151: PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1152: if (ni > 128) PetscCall(PetscFree(lix));
1153: } else {
1154: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1155: }
1156: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1157: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1158: PetscFunctionReturn(PETSC_SUCCESS);
1159: }
1161: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1162: {
1163: PetscFunctionBegin;
1166: VecCheckAssembled(x);
1168: if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1169: PetscAssertPointer(y, 3);
1170: for (PetscInt i = 0; i < nv; ++i) {
1173: PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1174: VecCheckSameSize(x, 1, y[i], 3);
1175: VecCheckAssembled(y[i]);
1176: PetscCall(VecLockReadPush(y[i]));
1177: }
1178: PetscAssertPointer(result, 4);
1181: PetscCall(VecLockReadPush(x));
1182: PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1183: PetscCall((*mxdot)(x, nv, y, result));
1184: PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1185: PetscCall(VecLockReadPop(x));
1186: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1187: PetscFunctionReturn(PETSC_SUCCESS);
1188: }
1190: /*@
1191: VecMTDot - Computes indefinite vector multiple dot products.
1192: That is, it does NOT use the complex conjugate.
1194: Collective
1196: Input Parameters:
1197: + x - one vector
1198: . nv - number of vectors
1199: - y - array of vectors. Note that vectors are pointers
1201: Output Parameter:
1202: . val - array of the dot products
1204: Level: intermediate
1206: Notes for Users of Complex Numbers:
1207: For complex vectors, `VecMTDot()` computes the indefinite form
1208: $ val = (x,y) = y^T x,
1209: where y^T denotes the transpose of y.
1211: Use `VecMDot()` for the inner product
1212: $ val = (x,y) = y^H x,
1213: where y^H denotes the conjugate transpose of y.
1215: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1216: @*/
1217: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1218: {
1219: PetscFunctionBegin;
1221: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1222: PetscFunctionReturn(PETSC_SUCCESS);
1223: }
1225: /*@
1226: VecMDot - Computes multiple vector dot products.
1228: Collective
1230: Input Parameters:
1231: + x - one vector
1232: . nv - number of vectors
1233: - y - array of vectors.
1235: Output Parameter:
1236: . val - array of the dot products (does not allocate the array)
1238: Level: intermediate
1240: Notes for Users of Complex Numbers:
1241: For complex vectors, `VecMDot()` computes
1242: $ val = (x,y) = y^H x,
1243: where y^H denotes the conjugate transpose of y.
1245: Use `VecMTDot()` for the indefinite form
1246: $ val = (x,y) = y^T x,
1247: where y^T denotes the transpose of y.
1249: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`
1250: @*/
1251: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1252: {
1253: PetscFunctionBegin;
1255: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1256: PetscFunctionReturn(PETSC_SUCCESS);
1257: }
1259: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
1260: {
1261: PetscFunctionBegin;
1263: VecCheckAssembled(y);
1265: PetscCall(VecSetErrorIfLocked(y, 1));
1266: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1267: if (nv) {
1268: PetscInt zeros = 0;
1270: PetscAssertPointer(alpha, 3);
1271: PetscAssertPointer(x, 4);
1272: for (PetscInt i = 0; i < nv; ++i) {
1276: PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1277: VecCheckSameSize(y, 1, x[i], 4);
1278: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1279: VecCheckAssembled(x[i]);
1280: PetscCall(VecLockReadPush(x[i]));
1281: zeros += alpha[i] == (PetscScalar)0.0;
1282: }
1284: if (zeros < nv) {
1285: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1286: VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1287: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1288: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1289: }
1291: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1292: }
1293: PetscFunctionReturn(PETSC_SUCCESS);
1294: }
1296: /*@
1297: VecMAXPY - Computes `y = y + sum alpha[i] x[i]`
1299: Logically Collective
1301: Input Parameters:
1302: + nv - number of scalars and x-vectors
1303: . alpha - array of scalars
1304: . y - one vector
1305: - x - array of vectors
1307: Level: intermediate
1309: Note:
1310: `y` cannot be any of the `x` vectors
1312: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`,`VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1313: @*/
1314: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1315: {
1316: PetscFunctionBegin;
1317: PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1318: PetscFunctionReturn(PETSC_SUCCESS);
1319: }
1321: /*@
1322: VecMAXPBY - Computes `y = beta y + sum alpha[i] x[i]`
1324: Logically Collective
1326: Input Parameters:
1327: + nv - number of scalars and x-vectors
1328: . alpha - array of scalars
1329: . beta - scalar
1330: . y - one vector
1331: - x - array of vectors
1333: Level: intermediate
1335: Note:
1336: `y` cannot be any of the `x` vectors.
1338: Developer Notes:
1339: This is a convenience routine, but implementations might be able to optimize it, for example, when `beta` is zero.
1341: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1342: @*/
1343: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
1344: {
1345: PetscFunctionBegin;
1347: VecCheckAssembled(y);
1349: PetscCall(VecSetErrorIfLocked(y, 1));
1350: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1353: if (y->ops->maxpby) {
1354: PetscInt zeros = 0;
1356: if (nv) {
1357: PetscAssertPointer(alpha, 3);
1358: PetscAssertPointer(x, 5);
1359: }
1361: for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
1365: PetscCheckSameTypeAndComm(y, 1, x[i], 5);
1366: VecCheckSameSize(y, 1, x[i], 5);
1367: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1368: VecCheckAssembled(x[i]);
1369: PetscCall(VecLockReadPush(x[i]));
1370: zeros += alpha[i] == (PetscScalar)0.0;
1371: }
1373: if (zeros < nv) { // has nonzero alpha
1374: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1375: PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1376: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1377: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1378: } else {
1379: PetscCall(VecScale(y, beta));
1380: }
1382: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1383: } else { // no maxpby
1384: if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1385: else PetscCall(VecScale(y, beta));
1386: PetscCall(VecMAXPY(y, nv, alpha, x));
1387: }
1388: PetscFunctionReturn(PETSC_SUCCESS);
1389: }
1391: /*@
1392: VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
1393: in the order they appear in the array. The concatenated vector resides on the same
1394: communicator and is the same type as the source vectors.
1396: Collective
1398: Input Parameters:
1399: + nx - number of vectors to be concatenated
1400: - X - array containing the vectors to be concatenated in the order of concatenation
1402: Output Parameters:
1403: + Y - concatenated vector
1404: - x_is - array of index sets corresponding to the concatenated components of `Y` (pass `NULL` if not needed)
1406: Level: advanced
1408: Notes:
1409: Concatenation is similar to the functionality of a `VECNEST` object; they both represent combination of
1410: different vector spaces. However, concatenated vectors do not store any information about their
1411: sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1412: manipulation of data in the concatenated vector that corresponds to the original components at creation.
1414: This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1415: has to operate on combined vector spaces and cannot utilize `VECNEST` objects due to incompatibility with
1416: bound projections.
1418: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1419: @*/
1420: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1421: {
1422: MPI_Comm comm;
1423: VecType vec_type;
1424: Vec Ytmp, Xtmp;
1425: IS *is_tmp;
1426: PetscInt i, shift = 0, Xnl, Xng, Xbegin;
1428: PetscFunctionBegin;
1432: PetscAssertPointer(Y, 3);
1434: if ((*X)->ops->concatenate) {
1435: /* use the dedicated concatenation function if available */
1436: PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1437: } else {
1438: /* loop over vectors and start creating IS */
1439: comm = PetscObjectComm((PetscObject)*X);
1440: PetscCall(VecGetType(*X, &vec_type));
1441: PetscCall(PetscMalloc1(nx, &is_tmp));
1442: for (i = 0; i < nx; i++) {
1443: PetscCall(VecGetSize(X[i], &Xng));
1444: PetscCall(VecGetLocalSize(X[i], &Xnl));
1445: PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1446: PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1447: shift += Xng;
1448: }
1449: /* create the concatenated vector */
1450: PetscCall(VecCreate(comm, &Ytmp));
1451: PetscCall(VecSetType(Ytmp, vec_type));
1452: PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1453: PetscCall(VecSetUp(Ytmp));
1454: /* copy data from X array to Y and return */
1455: for (i = 0; i < nx; i++) {
1456: PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1457: PetscCall(VecCopy(X[i], Xtmp));
1458: PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1459: }
1460: *Y = Ytmp;
1461: if (x_is) {
1462: *x_is = is_tmp;
1463: } else {
1464: for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1465: PetscCall(PetscFree(is_tmp));
1466: }
1467: }
1468: PetscFunctionReturn(PETSC_SUCCESS);
1469: }
1471: /* A helper function for VecGetSubVector to check if we can implement it with no-copy (i.e. the subvector shares
1472: memory with the original vector), and the block size of the subvector.
1474: Input Parameters:
1475: + X - the original vector
1476: - is - the index set of the subvector
1478: Output Parameters:
1479: + contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
1480: . start - start of contiguous block, as an offset from the start of the ownership range of the original vector
1481: - blocksize - the block size of the subvector
1483: */
1484: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1485: {
1486: PetscInt gstart, gend, lstart;
1487: PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1488: PetscInt n, N, ibs, vbs, bs = 1;
1490: PetscFunctionBegin;
1491: PetscCall(ISGetLocalSize(is, &n));
1492: PetscCall(ISGetSize(is, &N));
1493: PetscCall(ISGetBlockSize(is, &ibs));
1494: PetscCall(VecGetBlockSize(X, &vbs));
1495: PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1496: PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1497: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1498: if (ibs > 1) {
1499: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1500: bs = ibs;
1501: } else {
1502: if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1503: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1504: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1505: }
1507: *contig = red[0];
1508: *start = lstart;
1509: *blocksize = bs;
1510: PetscFunctionReturn(PETSC_SUCCESS);
1511: }
1513: /* A helper function for VecGetSubVector, to be used when we have to build a standalone subvector through VecScatter
1515: Input Parameters:
1516: + X - the original vector
1517: . is - the index set of the subvector
1518: - bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()
1520: Output Parameter:
1521: . Z - the subvector, which will compose the VecScatter context on output
1522: */
1523: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1524: {
1525: PetscInt n, N;
1526: VecScatter vscat;
1527: Vec Y;
1529: PetscFunctionBegin;
1530: PetscCall(ISGetLocalSize(is, &n));
1531: PetscCall(ISGetSize(is, &N));
1532: PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1533: PetscCall(VecSetSizes(Y, n, N));
1534: PetscCall(VecSetBlockSize(Y, bs));
1535: PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1536: PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1537: PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1538: PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1539: PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1540: PetscCall(VecScatterDestroy(&vscat));
1541: *Z = Y;
1542: PetscFunctionReturn(PETSC_SUCCESS);
1543: }
1545: /*@
1546: VecGetSubVector - Gets a vector representing part of another vector
1548: Collective
1550: Input Parameters:
1551: + X - vector from which to extract a subvector
1552: - is - index set representing portion of `X` to extract
1554: Output Parameter:
1555: . Y - subvector corresponding to `is`
1557: Level: advanced
1559: Notes:
1560: The subvector `Y` should be returned with `VecRestoreSubVector()`.
1561: `X` and `is` must be defined on the same communicator
1563: Changes to the subvector will be reflected in the `X` vector on the call to `VecRestoreSubVector()`.
1565: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1566: modifying the subvector. Other non-overlapping subvectors can still be obtained from `X` using this function.
1568: The resulting subvector inherits the block size from `is` if greater than one. Otherwise, the block size is guessed from the block size of the original `X`.
1570: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1571: @*/
1572: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1573: {
1574: Vec Z;
1576: PetscFunctionBegin;
1579: PetscCheckSameComm(X, 1, is, 2);
1580: PetscAssertPointer(Y, 3);
1581: if (X->ops->getsubvector) {
1582: PetscUseTypeMethod(X, getsubvector, is, &Z);
1583: } else { /* Default implementation currently does no caching */
1584: PetscBool contig;
1585: PetscInt n, N, start, bs;
1587: PetscCall(ISGetLocalSize(is, &n));
1588: PetscCall(ISGetSize(is, &N));
1589: PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1590: if (contig) { /* We can do a no-copy implementation */
1591: const PetscScalar *x;
1592: PetscInt state = 0;
1593: PetscBool isstd, iscuda, iship;
1595: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1596: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1597: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1598: if (iscuda) {
1599: #if defined(PETSC_HAVE_CUDA)
1600: const PetscScalar *x_d;
1601: PetscMPIInt size;
1602: PetscOffloadMask flg;
1604: PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1605: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1606: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1607: if (x) x += start;
1608: if (x_d) x_d += start;
1609: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1610: if (size == 1) {
1611: PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1612: } else {
1613: PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1614: }
1615: Z->offloadmask = flg;
1616: #endif
1617: } else if (iship) {
1618: #if defined(PETSC_HAVE_HIP)
1619: const PetscScalar *x_d;
1620: PetscMPIInt size;
1621: PetscOffloadMask flg;
1623: PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1624: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1625: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1626: if (x) x += start;
1627: if (x_d) x_d += start;
1628: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1629: if (size == 1) {
1630: PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1631: } else {
1632: PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1633: }
1634: Z->offloadmask = flg;
1635: #endif
1636: } else if (isstd) {
1637: PetscMPIInt size;
1639: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1640: PetscCall(VecGetArrayRead(X, &x));
1641: if (x) x += start;
1642: if (size == 1) {
1643: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1644: } else {
1645: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1646: }
1647: PetscCall(VecRestoreArrayRead(X, &x));
1648: } else { /* default implementation: use place array */
1649: PetscCall(VecGetArrayRead(X, &x));
1650: PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1651: PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1652: PetscCall(VecSetSizes(Z, n, N));
1653: PetscCall(VecSetBlockSize(Z, bs));
1654: PetscCall(VecPlaceArray(Z, PetscSafePointerPlusOffset(x, start)));
1655: PetscCall(VecRestoreArrayRead(X, &x));
1656: }
1658: /* this is relevant only in debug mode */
1659: PetscCall(VecLockGet(X, &state));
1660: if (state) PetscCall(VecLockReadPush(Z));
1661: Z->ops->placearray = NULL;
1662: Z->ops->replacearray = NULL;
1663: } else { /* Have to create a scatter and do a copy */
1664: PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1665: }
1666: }
1667: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1668: if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1669: PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1670: *Y = Z;
1671: PetscFunctionReturn(PETSC_SUCCESS);
1672: }
1674: /*@
1675: VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`
1677: Collective
1679: Input Parameters:
1680: + X - vector from which subvector was obtained
1681: . is - index set representing the subset of `X`
1682: - Y - subvector being restored
1684: Level: advanced
1686: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1687: @*/
1688: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1689: {
1690: PETSC_UNUSED PetscObjectState dummystate = 0;
1691: PetscBool unchanged;
1693: PetscFunctionBegin;
1696: PetscCheckSameComm(X, 1, is, 2);
1697: PetscAssertPointer(Y, 3);
1700: if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1701: else {
1702: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1703: if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1704: VecScatter scatter;
1705: PetscInt state;
1707: PetscCall(VecLockGet(X, &state));
1708: PetscCheck(state == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec X is locked for read-only or read/write access");
1710: PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1711: if (scatter) {
1712: PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1713: PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1714: } else {
1715: PetscBool iscuda, iship;
1716: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1717: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1719: if (iscuda) {
1720: #if defined(PETSC_HAVE_CUDA)
1721: PetscOffloadMask ymask = (*Y)->offloadmask;
1723: /* The offloadmask of X dictates where to move memory
1724: If X GPU data is valid, then move Y data on GPU if needed
1725: Otherwise, move back to the CPU */
1726: switch (X->offloadmask) {
1727: case PETSC_OFFLOAD_BOTH:
1728: if (ymask == PETSC_OFFLOAD_CPU) {
1729: PetscCall(VecCUDAResetArray(*Y));
1730: } else if (ymask == PETSC_OFFLOAD_GPU) {
1731: X->offloadmask = PETSC_OFFLOAD_GPU;
1732: }
1733: break;
1734: case PETSC_OFFLOAD_GPU:
1735: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1736: break;
1737: case PETSC_OFFLOAD_CPU:
1738: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1739: break;
1740: case PETSC_OFFLOAD_UNALLOCATED:
1741: case PETSC_OFFLOAD_KOKKOS:
1742: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1743: }
1744: #endif
1745: } else if (iship) {
1746: #if defined(PETSC_HAVE_HIP)
1747: PetscOffloadMask ymask = (*Y)->offloadmask;
1749: /* The offloadmask of X dictates where to move memory
1750: If X GPU data is valid, then move Y data on GPU if needed
1751: Otherwise, move back to the CPU */
1752: switch (X->offloadmask) {
1753: case PETSC_OFFLOAD_BOTH:
1754: if (ymask == PETSC_OFFLOAD_CPU) {
1755: PetscCall(VecHIPResetArray(*Y));
1756: } else if (ymask == PETSC_OFFLOAD_GPU) {
1757: X->offloadmask = PETSC_OFFLOAD_GPU;
1758: }
1759: break;
1760: case PETSC_OFFLOAD_GPU:
1761: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1762: break;
1763: case PETSC_OFFLOAD_CPU:
1764: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1765: break;
1766: case PETSC_OFFLOAD_UNALLOCATED:
1767: case PETSC_OFFLOAD_KOKKOS:
1768: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1769: }
1770: #endif
1771: } else {
1772: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1773: PetscCall(VecResetArray(*Y));
1774: }
1775: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1776: }
1777: }
1778: }
1779: PetscCall(VecDestroy(Y));
1780: PetscFunctionReturn(PETSC_SUCCESS);
1781: }
1783: /*@
1784: VecCreateLocalVector - Creates a vector object suitable for use with `VecGetLocalVector()` and friends. You must call `VecDestroy()` when the
1785: vector is no longer needed.
1787: Not Collective.
1789: Input Parameter:
1790: . v - The vector for which the local vector is desired.
1792: Output Parameter:
1793: . w - Upon exit this contains the local vector.
1795: Level: beginner
1797: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1798: @*/
1799: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1800: {
1801: VecType roottype;
1802: PetscInt n;
1804: PetscFunctionBegin;
1806: PetscAssertPointer(w, 2);
1807: if (v->ops->createlocalvector) {
1808: PetscUseTypeMethod(v, createlocalvector, w);
1809: PetscFunctionReturn(PETSC_SUCCESS);
1810: }
1811: PetscCall(VecGetRootType_Private(v, &roottype));
1812: PetscCall(VecCreate(PETSC_COMM_SELF, w));
1813: PetscCall(VecGetLocalSize(v, &n));
1814: PetscCall(VecSetSizes(*w, n, n));
1815: PetscCall(VecGetBlockSize(v, &n));
1816: PetscCall(VecSetBlockSize(*w, n));
1817: PetscCall(VecSetType(*w, roottype));
1818: PetscFunctionReturn(PETSC_SUCCESS);
1819: }
1821: /*@
1822: VecGetLocalVectorRead - Maps the local portion of a vector into a
1823: vector.
1825: Not Collective.
1827: Input Parameter:
1828: . v - The vector for which the local vector is desired.
1830: Output Parameter:
1831: . w - Upon exit this contains the local vector.
1833: Level: beginner
1835: Notes:
1836: You must call `VecRestoreLocalVectorRead()` when the local
1837: vector is no longer needed.
1839: This function is similar to `VecGetArrayRead()` which maps the local
1840: portion into a raw pointer. `VecGetLocalVectorRead()` is usually
1841: almost as efficient as `VecGetArrayRead()` but in certain circumstances
1842: `VecGetLocalVectorRead()` can be much more efficient than
1843: `VecGetArrayRead()`. This is because the construction of a contiguous
1844: array representing the vector data required by `VecGetArrayRead()` can
1845: be an expensive operation for certain vector types. For example, for
1846: GPU vectors `VecGetArrayRead()` requires that the data between device
1847: and host is synchronized.
1849: Unlike `VecGetLocalVector()`, this routine is not collective and
1850: preserves cached information.
1852: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1853: @*/
1854: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1855: {
1856: PetscFunctionBegin;
1859: VecCheckSameLocalSize(v, 1, w, 2);
1860: if (v->ops->getlocalvectorread) {
1861: PetscUseTypeMethod(v, getlocalvectorread, w);
1862: } else {
1863: PetscScalar *a;
1865: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1866: PetscCall(VecPlaceArray(w, a));
1867: }
1868: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1869: PetscCall(VecLockReadPush(v));
1870: PetscCall(VecLockReadPush(w));
1871: PetscFunctionReturn(PETSC_SUCCESS);
1872: }
1874: /*@
1875: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1876: previously mapped into a vector using `VecGetLocalVectorRead()`.
1878: Not Collective.
1880: Input Parameters:
1881: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
1882: - w - The vector into which the local portion of `v` was mapped.
1884: Level: beginner
1886: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1887: @*/
1888: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1889: {
1890: PetscFunctionBegin;
1893: if (v->ops->restorelocalvectorread) {
1894: PetscUseTypeMethod(v, restorelocalvectorread, w);
1895: } else {
1896: const PetscScalar *a;
1898: PetscCall(VecGetArrayRead(w, &a));
1899: PetscCall(VecRestoreArrayRead(v, &a));
1900: PetscCall(VecResetArray(w));
1901: }
1902: PetscCall(VecLockReadPop(v));
1903: PetscCall(VecLockReadPop(w));
1904: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1905: PetscFunctionReturn(PETSC_SUCCESS);
1906: }
1908: /*@
1909: VecGetLocalVector - Maps the local portion of a vector into a
1910: vector.
1912: Collective
1914: Input Parameter:
1915: . v - The vector for which the local vector is desired.
1917: Output Parameter:
1918: . w - Upon exit this contains the local vector.
1920: Level: beginner
1922: Notes:
1923: You must call `VecRestoreLocalVector()` when the local
1924: vector is no longer needed.
1926: This function is similar to `VecGetArray()` which maps the local
1927: portion into a raw pointer. `VecGetLocalVector()` is usually about as
1928: efficient as `VecGetArray()` but in certain circumstances
1929: `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1930: This is because the construction of a contiguous array representing
1931: the vector data required by `VecGetArray()` can be an expensive
1932: operation for certain vector types. For example, for GPU vectors
1933: `VecGetArray()` requires that the data between device and host is
1934: synchronized.
1936: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1937: @*/
1938: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1939: {
1940: PetscFunctionBegin;
1943: VecCheckSameLocalSize(v, 1, w, 2);
1944: if (v->ops->getlocalvector) {
1945: PetscUseTypeMethod(v, getlocalvector, w);
1946: } else {
1947: PetscScalar *a;
1949: PetscCall(VecGetArray(v, &a));
1950: PetscCall(VecPlaceArray(w, a));
1951: }
1952: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1953: PetscFunctionReturn(PETSC_SUCCESS);
1954: }
1956: /*@
1957: VecRestoreLocalVector - Unmaps the local portion of a vector
1958: previously mapped into a vector using `VecGetLocalVector()`.
1960: Logically Collective.
1962: Input Parameters:
1963: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
1964: - w - The vector into which the local portion of `v` was mapped.
1966: Level: beginner
1968: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1969: @*/
1970: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1971: {
1972: PetscFunctionBegin;
1975: if (v->ops->restorelocalvector) {
1976: PetscUseTypeMethod(v, restorelocalvector, w);
1977: } else {
1978: PetscScalar *a;
1979: PetscCall(VecGetArray(w, &a));
1980: PetscCall(VecRestoreArray(v, &a));
1981: PetscCall(VecResetArray(w));
1982: }
1983: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1984: PetscCall(PetscObjectStateIncrease((PetscObject)v));
1985: PetscFunctionReturn(PETSC_SUCCESS);
1986: }
1988: /*@C
1989: VecGetArray - Returns a pointer to a contiguous array that contains this
1990: MPI processes's portion of the vector data
1992: Logically Collective
1994: Input Parameter:
1995: . x - the vector
1997: Output Parameter:
1998: . a - location to put pointer to the array
2000: Level: beginner
2002: Notes:
2003: For the standard PETSc vectors, `VecGetArray()` returns a pointer to the local data array and
2004: does not use any copies. If the underlying vector data is not stored in a contiguous array
2005: this routine will copy the data to a contiguous array and return a pointer to that. You MUST
2006: call `VecRestoreArray()` when you no longer need access to the array.
2008: Fortran Note:
2009: .vb
2010: PetscScalar, pointer :: a(:)
2011: .ve
2013: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecPlaceArray()`, `VecGetArray2d()`,
2014: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2015: @*/
2016: PetscErrorCode VecGetArray(Vec x, PetscScalar *a[])
2017: {
2018: PetscFunctionBegin;
2020: PetscCall(VecSetErrorIfLocked(x, 1));
2021: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
2022: PetscUseTypeMethod(x, getarray, a);
2023: } else if (x->petscnative) { /* VECSTANDARD */
2024: *a = *((PetscScalar **)x->data);
2025: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
2026: PetscFunctionReturn(PETSC_SUCCESS);
2027: }
2029: /*@C
2030: VecRestoreArray - Restores a vector after `VecGetArray()` has been called and the array is no longer needed
2032: Logically Collective
2034: Input Parameters:
2035: + x - the vector
2036: - a - location of pointer to array obtained from `VecGetArray()`
2038: Level: beginner
2040: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2041: `VecGetArrayPair()`, `VecRestoreArrayPair()`
2042: @*/
2043: PetscErrorCode VecRestoreArray(Vec x, PetscScalar *a[])
2044: {
2045: PetscFunctionBegin;
2047: if (a) PetscAssertPointer(a, 2);
2048: if (x->ops->restorearray) {
2049: PetscUseTypeMethod(x, restorearray, a);
2050: } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
2051: if (a) *a = NULL;
2052: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2053: PetscFunctionReturn(PETSC_SUCCESS);
2054: }
2055: /*@C
2056: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
2058: Not Collective
2060: Input Parameter:
2061: . x - the vector
2063: Output Parameter:
2064: . a - the array
2066: Level: beginner
2068: Notes:
2069: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2071: Unlike `VecGetArray()`, preserves cached information like vector norms.
2073: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
2074: implementations may require a copy, but such implementations should cache the contiguous representation so that
2075: only one copy is performed when this routine is called multiple times in sequence.
2077: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2078: @*/
2079: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar *a[])
2080: {
2081: PetscFunctionBegin;
2083: PetscAssertPointer(a, 2);
2084: if (x->ops->getarrayread) {
2085: PetscUseTypeMethod(x, getarrayread, a);
2086: } else if (x->ops->getarray) {
2087: PetscObjectState state;
2089: /* VECNEST, VECCUDA, VECKOKKOS etc */
2090: // x->ops->getarray may bump the object state, but since we know this is a read-only get
2091: // we can just undo that
2092: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2093: PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
2094: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2095: } else if (x->petscnative) {
2096: /* VECSTANDARD */
2097: *a = *((PetscScalar **)x->data);
2098: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
2099: PetscFunctionReturn(PETSC_SUCCESS);
2100: }
2102: /*@C
2103: VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`
2105: Not Collective
2107: Input Parameters:
2108: + x - the vector
2109: - a - the array
2111: Level: beginner
2113: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2114: @*/
2115: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar *a[])
2116: {
2117: PetscFunctionBegin;
2119: if (a) PetscAssertPointer(a, 2);
2120: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
2121: /* nothing */
2122: } else if (x->ops->restorearrayread) { /* VECNEST */
2123: PetscUseTypeMethod(x, restorearrayread, a);
2124: } else { /* No one? */
2125: PetscObjectState state;
2127: // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2128: // we can just undo that
2129: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2130: PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2131: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2132: }
2133: if (a) *a = NULL;
2134: PetscFunctionReturn(PETSC_SUCCESS);
2135: }
2137: /*@C
2138: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
2139: MPI processes's portion of the vector data.
2141: Logically Collective
2143: Input Parameter:
2144: . x - the vector
2146: Output Parameter:
2147: . a - location to put pointer to the array
2149: Level: intermediate
2151: Note:
2152: The values in this array are NOT valid, the caller of this routine is responsible for putting
2153: values into the array; any values it does not set will be invalid.
2155: The array must be returned using a matching call to `VecRestoreArrayWrite()`.
2157: For vectors associated with GPUs, the host and device vectors are not synchronized before
2158: giving access. If you need correct values in the array use `VecGetArray()`
2160: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecPlaceArray()`, `VecGetArray2d()`,
2161: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
2162: @*/
2163: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar *a[])
2164: {
2165: PetscFunctionBegin;
2167: PetscAssertPointer(a, 2);
2168: PetscCall(VecSetErrorIfLocked(x, 1));
2169: if (x->ops->getarraywrite) {
2170: PetscUseTypeMethod(x, getarraywrite, a);
2171: } else {
2172: PetscCall(VecGetArray(x, a));
2173: }
2174: PetscFunctionReturn(PETSC_SUCCESS);
2175: }
2177: /*@C
2178: VecRestoreArrayWrite - Restores a vector after `VecGetArrayWrite()` has been called.
2180: Logically Collective
2182: Input Parameters:
2183: + x - the vector
2184: - a - location of pointer to array obtained from `VecGetArray()`
2186: Level: beginner
2188: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2189: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2190: @*/
2191: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar *a[])
2192: {
2193: PetscFunctionBegin;
2195: if (a) PetscAssertPointer(a, 2);
2196: if (x->ops->restorearraywrite) {
2197: PetscUseTypeMethod(x, restorearraywrite, a);
2198: } else if (x->ops->restorearray) {
2199: PetscUseTypeMethod(x, restorearray, a);
2200: }
2201: if (a) *a = NULL;
2202: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2203: PetscFunctionReturn(PETSC_SUCCESS);
2204: }
2206: /*@C
2207: VecGetArrays - Returns a pointer to the arrays in a set of vectors
2208: that were created by a call to `VecDuplicateVecs()`.
2210: Logically Collective; No Fortran Support
2212: Input Parameters:
2213: + x - the vectors
2214: - n - the number of vectors
2216: Output Parameter:
2217: . a - location to put pointer to the array
2219: Level: intermediate
2221: Note:
2222: You MUST call `VecRestoreArrays()` when you no longer need access to the arrays.
2224: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2225: @*/
2226: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2227: {
2228: PetscInt i;
2229: PetscScalar **q;
2231: PetscFunctionBegin;
2232: PetscAssertPointer(x, 1);
2234: PetscAssertPointer(a, 3);
2235: PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2236: PetscCall(PetscMalloc1(n, &q));
2237: for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2238: *a = q;
2239: PetscFunctionReturn(PETSC_SUCCESS);
2240: }
2242: /*@C
2243: VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2244: has been called.
2246: Logically Collective; No Fortran Support
2248: Input Parameters:
2249: + x - the vector
2250: . n - the number of vectors
2251: - a - location of pointer to arrays obtained from `VecGetArrays()`
2253: Notes:
2254: For regular PETSc vectors this routine does not involve any copies. For
2255: any special vectors that do not store local vector data in a contiguous
2256: array, this routine will copy the data back into the underlying
2257: vector data structure from the arrays obtained with `VecGetArrays()`.
2259: Level: intermediate
2261: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2262: @*/
2263: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2264: {
2265: PetscInt i;
2266: PetscScalar **q = *a;
2268: PetscFunctionBegin;
2269: PetscAssertPointer(x, 1);
2271: PetscAssertPointer(a, 3);
2273: for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2274: PetscCall(PetscFree(q));
2275: PetscFunctionReturn(PETSC_SUCCESS);
2276: }
2278: /*@C
2279: VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g.,
2280: `VECCUDA`), the returned pointer will be a device pointer to the device memory that contains
2281: this MPI processes's portion of the vector data.
2283: Logically Collective; No Fortran Support
2285: Input Parameter:
2286: . x - the vector
2288: Output Parameters:
2289: + a - location to put pointer to the array
2290: - mtype - memory type of the array
2292: Level: beginner
2294: Note:
2295: Device data is guaranteed to have the latest value. Otherwise, when this is a host vector
2296: (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host
2297: pointer.
2299: For `VECKOKKOS`, if Kokkos is configured without device (e.g., use serial or openmp), per
2300: this function, the vector works like `VECSEQ`/`VECMPI`; otherwise, it works like `VECCUDA` or
2301: `VECHIP` etc.
2303: Use `VecRestoreArrayAndMemType()` when the array access is no longer needed.
2305: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`,
2306: `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2307: @*/
2308: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar *a[], PetscMemType *mtype)
2309: {
2310: PetscFunctionBegin;
2313: PetscAssertPointer(a, 2);
2314: if (mtype) PetscAssertPointer(mtype, 3);
2315: PetscCall(VecSetErrorIfLocked(x, 1));
2316: if (x->ops->getarrayandmemtype) {
2317: /* VECCUDA, VECKOKKOS etc */
2318: PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2319: } else {
2320: /* VECSTANDARD, VECNEST, VECVIENNACL */
2321: PetscCall(VecGetArray(x, a));
2322: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2323: }
2324: PetscFunctionReturn(PETSC_SUCCESS);
2325: }
2327: /*@C
2328: VecRestoreArrayAndMemType - Restores a vector after `VecGetArrayAndMemType()` has been called.
2330: Logically Collective; No Fortran Support
2332: Input Parameters:
2333: + x - the vector
2334: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`
2336: Level: beginner
2338: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`,
2339: `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2340: @*/
2341: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar *a[])
2342: {
2343: PetscFunctionBegin;
2346: if (a) PetscAssertPointer(a, 2);
2347: if (x->ops->restorearrayandmemtype) {
2348: /* VECCUDA, VECKOKKOS etc */
2349: PetscUseTypeMethod(x, restorearrayandmemtype, a);
2350: } else {
2351: /* VECNEST, VECVIENNACL */
2352: PetscCall(VecRestoreArray(x, a));
2353: } /* VECSTANDARD does nothing */
2354: if (a) *a = NULL;
2355: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2356: PetscFunctionReturn(PETSC_SUCCESS);
2357: }
2359: /*@C
2360: VecGetArrayReadAndMemType - Like `VecGetArrayRead()`, but if the input vector is a device vector, it will return a read-only device pointer.
2361: The returned pointer is guaranteed to point to up-to-date data. For host vectors, it functions as `VecGetArrayRead()`.
2363: Not Collective; No Fortran Support
2365: Input Parameter:
2366: . x - the vector
2368: Output Parameters:
2369: + a - the array
2370: - mtype - memory type of the array
2372: Level: beginner
2374: Notes:
2375: The array must be returned using a matching call to `VecRestoreArrayReadAndMemType()`.
2377: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2378: @*/
2379: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar *a[], PetscMemType *mtype)
2380: {
2381: PetscFunctionBegin;
2384: PetscAssertPointer(a, 2);
2385: if (mtype) PetscAssertPointer(mtype, 3);
2386: if (x->ops->getarrayreadandmemtype) {
2387: /* VECCUDA/VECHIP though they are also petscnative */
2388: PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2389: } else if (x->ops->getarrayandmemtype) {
2390: /* VECKOKKOS */
2391: PetscObjectState state;
2393: // see VecGetArrayRead() for why
2394: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2395: PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2396: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2397: } else {
2398: PetscCall(VecGetArrayRead(x, a));
2399: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2400: }
2401: PetscFunctionReturn(PETSC_SUCCESS);
2402: }
2404: /*@C
2405: VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`
2407: Not Collective; No Fortran Support
2409: Input Parameters:
2410: + x - the vector
2411: - a - the array
2413: Level: beginner
2415: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2416: @*/
2417: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar *a[])
2418: {
2419: PetscFunctionBegin;
2422: if (a) PetscAssertPointer(a, 2);
2423: if (x->ops->restorearrayreadandmemtype) {
2424: /* VECCUDA/VECHIP */
2425: PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2426: } else if (!x->petscnative) {
2427: /* VECNEST */
2428: PetscCall(VecRestoreArrayRead(x, a));
2429: }
2430: if (a) *a = NULL;
2431: PetscFunctionReturn(PETSC_SUCCESS);
2432: }
2434: /*@C
2435: VecGetArrayWriteAndMemType - Like `VecGetArrayWrite()`, but if this is a device vector it will always return
2436: a device pointer to the device memory that contains this processor's portion of the vector data.
2438: Logically Collective; No Fortran Support
2440: Input Parameter:
2441: . x - the vector
2443: Output Parameters:
2444: + a - the array
2445: - mtype - memory type of the array
2447: Level: beginner
2449: Note:
2450: The array must be returned using a matching call to `VecRestoreArrayWriteAndMemType()`, where it will label the device memory as most recent.
2452: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2453: @*/
2454: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar *a[], PetscMemType *mtype)
2455: {
2456: PetscFunctionBegin;
2459: PetscCall(VecSetErrorIfLocked(x, 1));
2460: PetscAssertPointer(a, 2);
2461: if (mtype) PetscAssertPointer(mtype, 3);
2462: if (x->ops->getarraywriteandmemtype) {
2463: /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2464: PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2465: } else if (x->ops->getarrayandmemtype) {
2466: PetscCall(VecGetArrayAndMemType(x, a, mtype));
2467: } else {
2468: /* VECNEST, VECVIENNACL */
2469: PetscCall(VecGetArrayWrite(x, a));
2470: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2471: }
2472: PetscFunctionReturn(PETSC_SUCCESS);
2473: }
2475: /*@C
2476: VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`
2478: Logically Collective; No Fortran Support
2480: Input Parameters:
2481: + x - the vector
2482: - a - the array
2484: Level: beginner
2486: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2487: @*/
2488: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar *a[])
2489: {
2490: PetscFunctionBegin;
2493: PetscCall(VecSetErrorIfLocked(x, 1));
2494: if (a) PetscAssertPointer(a, 2);
2495: if (x->ops->restorearraywriteandmemtype) {
2496: /* VECCUDA/VECHIP */
2497: PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2498: PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2499: } else if (x->ops->restorearrayandmemtype) {
2500: PetscCall(VecRestoreArrayAndMemType(x, a));
2501: } else {
2502: PetscCall(VecRestoreArray(x, a));
2503: }
2504: if (a) *a = NULL;
2505: PetscFunctionReturn(PETSC_SUCCESS);
2506: }
2508: /*@
2509: VecPlaceArray - Allows one to replace the array in a vector with an
2510: array provided by the user. This is useful to avoid copying an array
2511: into a vector.
2513: Logically Collective; No Fortran Support
2515: Input Parameters:
2516: + vec - the vector
2517: - array - the array
2519: Level: developer
2521: Notes:
2522: Adding `const` to `array` was an oversight, as subsequent operations on `vec` would
2523: likely modify the data in `array`. However, we have kept it to avoid breaking APIs.
2525: Use `VecReplaceArray()` instead to permanently replace the array
2527: You can return to the original array with a call to `VecResetArray()`. `vec` does not take
2528: ownership of `array` in any way.
2530: The user must free `array` themselves but be careful not to
2531: do so before the vector has either been destroyed, had its original array restored with
2532: `VecResetArray()` or permanently replaced with `VecReplaceArray()`.
2534: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2535: @*/
2536: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2537: {
2538: PetscFunctionBegin;
2541: if (array) PetscAssertPointer(array, 2);
2542: PetscUseTypeMethod(vec, placearray, array);
2543: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2544: PetscFunctionReturn(PETSC_SUCCESS);
2545: }
2547: /*@C
2548: VecReplaceArray - Allows one to replace the array in a vector with an
2549: array provided by the user. This is useful to avoid copying an array
2550: into a vector.
2552: Logically Collective; No Fortran Support
2554: Input Parameters:
2555: + vec - the vector
2556: - array - the array
2558: Level: developer
2560: Notes:
2561: Adding `const` to `array` was an oversight, as subsequent operations on `vec` would
2562: likely modify the data in `array`. However, we have kept it to avoid breaking APIs.
2564: This permanently replaces the array and frees the memory associated
2565: with the old array. Use `VecPlaceArray()` to temporarily replace the array.
2567: The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
2568: freed by the user. It will be freed when the vector is destroyed.
2570: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2571: @*/
2572: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2573: {
2574: PetscFunctionBegin;
2577: PetscUseTypeMethod(vec, replacearray, array);
2578: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2579: PetscFunctionReturn(PETSC_SUCCESS);
2580: }
2582: /*@C
2583: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2584: processor's portion of the vector data. You MUST call `VecRestoreArray2d()`
2585: when you no longer need access to the array.
2587: Logically Collective
2589: Input Parameters:
2590: + x - the vector
2591: . m - first dimension of two dimensional array
2592: . n - second dimension of 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)
2596: Output Parameter:
2597: . a - location to put pointer to the array
2599: Level: developer
2601: Notes:
2602: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2603: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2604: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2605: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2607: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2609: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2610: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2611: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2612: @*/
2613: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2614: {
2615: PetscInt i, N;
2616: PetscScalar *aa;
2618: PetscFunctionBegin;
2620: PetscAssertPointer(a, 6);
2622: PetscCall(VecGetLocalSize(x, &N));
2623: PetscCheck(m * n == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 2d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n);
2624: PetscCall(VecGetArray(x, &aa));
2626: PetscCall(PetscMalloc1(m, a));
2627: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2628: *a -= mstart;
2629: PetscFunctionReturn(PETSC_SUCCESS);
2630: }
2632: /*@C
2633: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2634: processor's portion of the vector data. You MUST call `VecRestoreArray2dWrite()`
2635: when you no longer need access to the array.
2637: Logically Collective
2639: Input Parameters:
2640: + x - the vector
2641: . m - first dimension of two dimensional array
2642: . n - second dimension of two dimensional array
2643: . mstart - first index you will use in first coordinate direction (often 0)
2644: - nstart - first index in the second coordinate direction (often 0)
2646: Output Parameter:
2647: . a - location to put pointer to the array
2649: Level: developer
2651: Notes:
2652: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2653: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2654: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2655: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2657: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2659: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2660: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2661: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2662: @*/
2663: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2664: {
2665: PetscInt i, N;
2666: PetscScalar *aa;
2668: PetscFunctionBegin;
2670: PetscAssertPointer(a, 6);
2672: PetscCall(VecGetLocalSize(x, &N));
2673: PetscCheck(m * n == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 2d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n);
2674: PetscCall(VecGetArrayWrite(x, &aa));
2676: PetscCall(PetscMalloc1(m, a));
2677: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2678: *a -= mstart;
2679: PetscFunctionReturn(PETSC_SUCCESS);
2680: }
2682: /*@C
2683: VecRestoreArray2d - Restores a vector after `VecGetArray2d()` has been called.
2685: Logically Collective
2687: Input Parameters:
2688: + x - the vector
2689: . m - first dimension of two dimensional array
2690: . n - second dimension of the two dimensional array
2691: . mstart - first index you will use in first coordinate direction (often 0)
2692: . nstart - first index in the second coordinate direction (often 0)
2693: - a - location of pointer to array obtained from `VecGetArray2d()`
2695: Level: developer
2697: Notes:
2698: For regular PETSc vectors this routine does not involve any copies. For
2699: any special vectors that do not store local vector data in a contiguous
2700: array, this routine will copy the data back into the underlying
2701: vector data structure from the array obtained with `VecGetArray()`.
2703: This routine actually zeros out the `a` pointer.
2705: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2706: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2707: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2708: @*/
2709: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2710: {
2711: void *dummy;
2713: PetscFunctionBegin;
2715: PetscAssertPointer(a, 6);
2717: dummy = (void *)(*a + mstart);
2718: PetscCall(PetscFree(dummy));
2719: PetscCall(VecRestoreArray(x, NULL));
2720: *a = NULL;
2721: PetscFunctionReturn(PETSC_SUCCESS);
2722: }
2724: /*@C
2725: VecRestoreArray2dWrite - Restores a vector after `VecGetArray2dWrite()` has been called.
2727: Logically Collective
2729: Input Parameters:
2730: + x - the vector
2731: . m - first dimension of two dimensional array
2732: . n - second dimension of the two dimensional array
2733: . mstart - first index you will use in first coordinate direction (often 0)
2734: . nstart - first index in the second coordinate direction (often 0)
2735: - a - location of pointer to array obtained from `VecGetArray2d()`
2737: Level: developer
2739: Notes:
2740: For regular PETSc vectors this routine does not involve any copies. For
2741: any special vectors that do not store local vector data in a contiguous
2742: array, this routine will copy the data back into the underlying
2743: vector data structure from the array obtained with `VecGetArray()`.
2745: This routine actually zeros out the `a` pointer.
2747: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2748: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2749: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2750: @*/
2751: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2752: {
2753: void *dummy;
2755: PetscFunctionBegin;
2757: PetscAssertPointer(a, 6);
2759: dummy = (void *)(*a + mstart);
2760: PetscCall(PetscFree(dummy));
2761: PetscCall(VecRestoreArrayWrite(x, NULL));
2762: PetscFunctionReturn(PETSC_SUCCESS);
2763: }
2765: /*@C
2766: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2767: processor's portion of the vector data. You MUST call `VecRestoreArray1d()`
2768: when you no longer need access to the array.
2770: Logically Collective
2772: Input Parameters:
2773: + x - the vector
2774: . m - first dimension of two dimensional array
2775: - mstart - first index you will use in first coordinate direction (often 0)
2777: Output Parameter:
2778: . a - location to put pointer to the array
2780: Level: developer
2782: Notes:
2783: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2784: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2785: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2787: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2789: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2790: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2791: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2792: @*/
2793: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2794: {
2795: PetscInt N;
2797: PetscFunctionBegin;
2799: PetscAssertPointer(a, 4);
2801: PetscCall(VecGetLocalSize(x, &N));
2802: PetscCheck(m == N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local array size %" PetscInt_FMT " does not match 1d array dimensions %" PetscInt_FMT, N, m);
2803: PetscCall(VecGetArray(x, a));
2804: *a -= mstart;
2805: PetscFunctionReturn(PETSC_SUCCESS);
2806: }
2808: /*@C
2809: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
2810: processor's portion of the vector data. You MUST call `VecRestoreArray1dWrite()`
2811: when you no longer need access to the array.
2813: Logically Collective
2815: Input Parameters:
2816: + x - the vector
2817: . m - first dimension of two dimensional array
2818: - mstart - first index you will use in first coordinate direction (often 0)
2820: Output Parameter:
2821: . a - location to put pointer to the array
2823: Level: developer
2825: Notes:
2826: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2827: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2828: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2830: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2832: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2833: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2834: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2835: @*/
2836: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2837: {
2838: PetscInt N;
2840: PetscFunctionBegin;
2842: PetscAssertPointer(a, 4);
2844: PetscCall(VecGetLocalSize(x, &N));
2845: PetscCheck(m == N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local array size %" PetscInt_FMT " does not match 1d array dimensions %" PetscInt_FMT, N, m);
2846: PetscCall(VecGetArrayWrite(x, a));
2847: *a -= mstart;
2848: PetscFunctionReturn(PETSC_SUCCESS);
2849: }
2851: /*@C
2852: VecRestoreArray1d - Restores a vector after `VecGetArray1d()` has been called.
2854: Logically Collective
2856: Input Parameters:
2857: + x - the vector
2858: . m - first dimension of two dimensional array
2859: . mstart - first index you will use in first coordinate direction (often 0)
2860: - a - location of pointer to array obtained from `VecGetArray1d()`
2862: Level: developer
2864: Notes:
2865: For regular PETSc vectors this routine does not involve any copies. For
2866: any special vectors that do not store local vector data in a contiguous
2867: array, this routine will copy the data back into the underlying
2868: vector data structure from the array obtained with `VecGetArray1d()`.
2870: This routine actually zeros out the `a` pointer.
2872: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2873: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2874: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2875: @*/
2876: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2877: {
2878: PetscFunctionBegin;
2881: PetscCall(VecRestoreArray(x, NULL));
2882: *a = NULL;
2883: PetscFunctionReturn(PETSC_SUCCESS);
2884: }
2886: /*@C
2887: VecRestoreArray1dWrite - Restores a vector after `VecGetArray1dWrite()` has been called.
2889: Logically Collective
2891: Input Parameters:
2892: + x - the vector
2893: . m - first dimension of two dimensional array
2894: . mstart - first index you will use in first coordinate direction (often 0)
2895: - a - location of pointer to array obtained from `VecGetArray1d()`
2897: Level: developer
2899: Notes:
2900: For regular PETSc vectors this routine does not involve any copies. For
2901: any special vectors that do not store local vector data in a contiguous
2902: array, this routine will copy the data back into the underlying
2903: vector data structure from the array obtained with `VecGetArray1d()`.
2905: This routine actually zeros out the `a` pointer.
2907: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2908: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2909: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2910: @*/
2911: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2912: {
2913: PetscFunctionBegin;
2916: PetscCall(VecRestoreArrayWrite(x, NULL));
2917: *a = NULL;
2918: PetscFunctionReturn(PETSC_SUCCESS);
2919: }
2921: /*@C
2922: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
2923: processor's portion of the vector data. You MUST call `VecRestoreArray3d()`
2924: when you no longer need access to the array.
2926: Logically Collective
2928: Input Parameters:
2929: + x - the vector
2930: . m - first dimension of three dimensional array
2931: . n - second dimension of three dimensional array
2932: . p - third dimension of three dimensional array
2933: . mstart - first index you will use in first coordinate direction (often 0)
2934: . nstart - first index in the second coordinate direction (often 0)
2935: - pstart - first index in the third coordinate direction (often 0)
2937: Output Parameter:
2938: . a - location to put pointer to the array
2940: Level: developer
2942: Notes:
2943: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
2944: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2945: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2946: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
2948: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2950: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2951: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
2952: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2953: @*/
2954: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
2955: {
2956: PetscInt i, N, j;
2957: PetscScalar *aa, **b;
2959: PetscFunctionBegin;
2961: PetscAssertPointer(a, 8);
2963: PetscCall(VecGetLocalSize(x, &N));
2964: PetscCheck(m * n * p == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 3d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n, p);
2965: PetscCall(VecGetArray(x, &aa));
2967: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
2968: b = (PetscScalar **)((*a) + m);
2969: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
2970: for (i = 0; i < m; i++)
2971: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset(aa, i * n * p + j * p - pstart);
2972: *a -= mstart;
2973: PetscFunctionReturn(PETSC_SUCCESS);
2974: }
2976: /*@C
2977: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
2978: processor's portion of the vector data. You MUST call `VecRestoreArray3dWrite()`
2979: when you no longer need access to the array.
2981: Logically Collective
2983: Input Parameters:
2984: + x - the vector
2985: . m - first dimension of three dimensional array
2986: . n - second dimension of three dimensional array
2987: . p - third dimension of three dimensional array
2988: . mstart - first index you will use in first coordinate direction (often 0)
2989: . nstart - first index in the second coordinate direction (often 0)
2990: - pstart - first index in the third coordinate direction (often 0)
2992: Output Parameter:
2993: . a - location to put pointer to the array
2995: Level: developer
2997: Notes:
2998: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
2999: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3000: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3001: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3003: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3005: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3006: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3007: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3008: @*/
3009: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3010: {
3011: PetscInt i, N, j;
3012: PetscScalar *aa, **b;
3014: PetscFunctionBegin;
3016: PetscAssertPointer(a, 8);
3018: PetscCall(VecGetLocalSize(x, &N));
3019: PetscCheck(m * n * p == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 3d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n, p);
3020: PetscCall(VecGetArrayWrite(x, &aa));
3022: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3023: b = (PetscScalar **)((*a) + m);
3024: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3025: for (i = 0; i < m; i++)
3026: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3028: *a -= mstart;
3029: PetscFunctionReturn(PETSC_SUCCESS);
3030: }
3032: /*@C
3033: VecRestoreArray3d - Restores a vector after `VecGetArray3d()` has been called.
3035: Logically Collective
3037: Input Parameters:
3038: + x - the vector
3039: . m - first dimension of three dimensional array
3040: . n - second dimension of the three dimensional array
3041: . p - third dimension of the three dimensional array
3042: . mstart - first index you will use in first coordinate direction (often 0)
3043: . nstart - first index in the second coordinate direction (often 0)
3044: . pstart - first index in the third coordinate direction (often 0)
3045: - a - location of pointer to array obtained from VecGetArray3d()
3047: Level: developer
3049: Notes:
3050: For regular PETSc vectors this routine does not involve any copies. For
3051: any special vectors that do not store local vector data in a contiguous
3052: array, this routine will copy the data back into the underlying
3053: vector data structure from the array obtained with `VecGetArray()`.
3055: This routine actually zeros out the `a` pointer.
3057: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3058: `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3059: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3060: @*/
3061: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3062: {
3063: void *dummy;
3065: PetscFunctionBegin;
3067: PetscAssertPointer(a, 8);
3069: dummy = (void *)(*a + mstart);
3070: PetscCall(PetscFree(dummy));
3071: PetscCall(VecRestoreArray(x, NULL));
3072: *a = NULL;
3073: PetscFunctionReturn(PETSC_SUCCESS);
3074: }
3076: /*@C
3077: VecRestoreArray3dWrite - Restores a vector after `VecGetArray3dWrite()` has been called.
3079: Logically Collective
3081: Input Parameters:
3082: + x - the vector
3083: . m - first dimension of three dimensional array
3084: . n - second dimension of the three dimensional array
3085: . p - third dimension of the three dimensional array
3086: . mstart - first index you will use in first coordinate direction (often 0)
3087: . nstart - first index in the second coordinate direction (often 0)
3088: . pstart - first index in the third coordinate direction (often 0)
3089: - a - location of pointer to array obtained from VecGetArray3d()
3091: Level: developer
3093: Notes:
3094: For regular PETSc vectors this routine does not involve any copies. For
3095: any special vectors that do not store local vector data in a contiguous
3096: array, this routine will copy the data back into the underlying
3097: vector data structure from the array obtained with `VecGetArray()`.
3099: This routine actually zeros out the `a` pointer.
3101: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3102: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3103: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3104: @*/
3105: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3106: {
3107: void *dummy;
3109: PetscFunctionBegin;
3111: PetscAssertPointer(a, 8);
3113: dummy = (void *)(*a + mstart);
3114: PetscCall(PetscFree(dummy));
3115: PetscCall(VecRestoreArrayWrite(x, NULL));
3116: *a = NULL;
3117: PetscFunctionReturn(PETSC_SUCCESS);
3118: }
3120: /*@C
3121: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3122: processor's portion of the vector data. You MUST call `VecRestoreArray4d()`
3123: when you no longer need access to the array.
3125: Logically Collective
3127: Input Parameters:
3128: + x - the vector
3129: . m - first dimension of four dimensional array
3130: . n - second dimension of four dimensional array
3131: . p - third dimension of four dimensional array
3132: . q - fourth dimension of four dimensional array
3133: . mstart - first index you will use in first coordinate direction (often 0)
3134: . nstart - first index in the second coordinate direction (often 0)
3135: . pstart - first index in the third coordinate direction (often 0)
3136: - qstart - first index in the fourth coordinate direction (often 0)
3138: Output Parameter:
3139: . a - location to put pointer to the array
3141: Level: developer
3143: Notes:
3144: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3145: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3146: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3147: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3149: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3151: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3152: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3153: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3154: @*/
3155: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3156: {
3157: PetscInt i, N, j, k;
3158: PetscScalar *aa, ***b, **c;
3160: PetscFunctionBegin;
3162: PetscAssertPointer(a, 10);
3164: PetscCall(VecGetLocalSize(x, &N));
3165: PetscCheck(m * n * p * q == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 4d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n, p, q);
3166: PetscCall(VecGetArray(x, &aa));
3168: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3169: b = (PetscScalar ***)((*a) + m);
3170: c = (PetscScalar **)(b + m * n);
3171: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3172: for (i = 0; i < m; i++)
3173: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3174: for (i = 0; i < m; i++)
3175: for (j = 0; j < n; j++)
3176: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3177: *a -= mstart;
3178: PetscFunctionReturn(PETSC_SUCCESS);
3179: }
3181: /*@C
3182: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3183: processor's portion of the vector data. You MUST call `VecRestoreArray4dWrite()`
3184: when you no longer need access to the array.
3186: Logically Collective
3188: Input Parameters:
3189: + x - the vector
3190: . m - first dimension of four dimensional array
3191: . n - second dimension of four dimensional array
3192: . p - third dimension of four dimensional array
3193: . q - fourth dimension of four dimensional array
3194: . mstart - first index you will use in first coordinate direction (often 0)
3195: . nstart - first index in the second coordinate direction (often 0)
3196: . pstart - first index in the third coordinate direction (often 0)
3197: - qstart - first index in the fourth coordinate direction (often 0)
3199: Output Parameter:
3200: . a - location to put pointer to the array
3202: Level: developer
3204: Notes:
3205: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3206: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3207: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3208: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3210: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3212: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3213: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3214: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3215: @*/
3216: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3217: {
3218: PetscInt i, N, j, k;
3219: PetscScalar *aa, ***b, **c;
3221: PetscFunctionBegin;
3223: PetscAssertPointer(a, 10);
3225: PetscCall(VecGetLocalSize(x, &N));
3226: PetscCheck(m * n * p * q == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 4d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n, p, q);
3227: PetscCall(VecGetArrayWrite(x, &aa));
3229: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3230: b = (PetscScalar ***)((*a) + m);
3231: c = (PetscScalar **)(b + m * n);
3232: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3233: for (i = 0; i < m; i++)
3234: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3235: for (i = 0; i < m; i++)
3236: for (j = 0; j < n; j++)
3237: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3238: *a -= mstart;
3239: PetscFunctionReturn(PETSC_SUCCESS);
3240: }
3242: /*@C
3243: VecRestoreArray4d - Restores a vector after `VecGetArray4d()` has been called.
3245: Logically Collective
3247: Input Parameters:
3248: + x - the vector
3249: . m - first dimension of four dimensional array
3250: . n - second dimension of the four dimensional array
3251: . p - third dimension of the four dimensional array
3252: . q - fourth dimension of the four dimensional array
3253: . mstart - first index you will use in first coordinate direction (often 0)
3254: . nstart - first index in the second coordinate direction (often 0)
3255: . pstart - first index in the third coordinate direction (often 0)
3256: . qstart - first index in the fourth coordinate direction (often 0)
3257: - a - location of pointer to array obtained from VecGetArray4d()
3259: Level: developer
3261: Notes:
3262: For regular PETSc vectors this routine does not involve any copies. For
3263: any special vectors that do not store local vector data in a contiguous
3264: array, this routine will copy the data back into the underlying
3265: vector data structure from the array obtained with `VecGetArray()`.
3267: This routine actually zeros out the `a` pointer.
3269: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3270: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3271: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`
3272: @*/
3273: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3274: {
3275: void *dummy;
3277: PetscFunctionBegin;
3279: PetscAssertPointer(a, 10);
3281: dummy = (void *)(*a + mstart);
3282: PetscCall(PetscFree(dummy));
3283: PetscCall(VecRestoreArray(x, NULL));
3284: *a = NULL;
3285: PetscFunctionReturn(PETSC_SUCCESS);
3286: }
3288: /*@C
3289: VecRestoreArray4dWrite - Restores a vector after `VecGetArray4dWrite()` has been called.
3291: Logically Collective
3293: Input Parameters:
3294: + x - the vector
3295: . m - first dimension of four dimensional array
3296: . n - second dimension of the four dimensional array
3297: . p - third dimension of the four dimensional array
3298: . q - fourth dimension of the four dimensional array
3299: . mstart - first index you will use in first coordinate direction (often 0)
3300: . nstart - first index in the second coordinate direction (often 0)
3301: . pstart - first index in the third coordinate direction (often 0)
3302: . qstart - first index in the fourth coordinate direction (often 0)
3303: - a - location of pointer to array obtained from `VecGetArray4d()`
3305: Level: developer
3307: Notes:
3308: For regular PETSc vectors this routine does not involve any copies. For
3309: any special vectors that do not store local vector data in a contiguous
3310: array, this routine will copy the data back into the underlying
3311: vector data structure from the array obtained with `VecGetArray()`.
3313: This routine actually zeros out the `a` pointer.
3315: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3316: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3317: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3318: @*/
3319: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3320: {
3321: void *dummy;
3323: PetscFunctionBegin;
3325: PetscAssertPointer(a, 10);
3327: dummy = (void *)(*a + mstart);
3328: PetscCall(PetscFree(dummy));
3329: PetscCall(VecRestoreArrayWrite(x, NULL));
3330: *a = NULL;
3331: PetscFunctionReturn(PETSC_SUCCESS);
3332: }
3334: /*@C
3335: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3336: processor's portion of the vector data. You MUST call `VecRestoreArray2dRead()`
3337: when you no longer need access to the array.
3339: Logically Collective
3341: Input Parameters:
3342: + x - the vector
3343: . m - first dimension of two dimensional array
3344: . n - second dimension of two dimensional array
3345: . mstart - first index you will use in first coordinate direction (often 0)
3346: - nstart - first index in the second coordinate direction (often 0)
3348: Output Parameter:
3349: . a - location to put pointer to the array
3351: Level: developer
3353: Notes:
3354: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
3355: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3356: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3357: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
3359: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3361: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3362: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3363: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3364: @*/
3365: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3366: {
3367: PetscInt i, N;
3368: const PetscScalar *aa;
3370: PetscFunctionBegin;
3372: PetscAssertPointer(a, 6);
3374: PetscCall(VecGetLocalSize(x, &N));
3375: PetscCheck(m * n == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 2d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n);
3376: PetscCall(VecGetArrayRead(x, &aa));
3378: PetscCall(PetscMalloc1(m, a));
3379: for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3380: *a -= mstart;
3381: PetscFunctionReturn(PETSC_SUCCESS);
3382: }
3384: /*@C
3385: VecRestoreArray2dRead - Restores a vector after `VecGetArray2dRead()` has been called.
3387: Logically Collective
3389: Input Parameters:
3390: + x - the vector
3391: . m - first dimension of two dimensional array
3392: . n - second dimension of the two dimensional array
3393: . mstart - first index you will use in first coordinate direction (often 0)
3394: . nstart - first index in the second coordinate direction (often 0)
3395: - a - location of pointer to array obtained from VecGetArray2d()
3397: Level: developer
3399: Notes:
3400: For regular PETSc vectors this routine does not involve any copies. For
3401: any special vectors that do not store local vector data in a contiguous
3402: array, this routine will copy the data back into the underlying
3403: vector data structure from the array obtained with `VecGetArray()`.
3405: This routine actually zeros out the `a` pointer.
3407: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3408: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3409: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3410: @*/
3411: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3412: {
3413: void *dummy;
3415: PetscFunctionBegin;
3417: PetscAssertPointer(a, 6);
3419: dummy = (void *)(*a + mstart);
3420: PetscCall(PetscFree(dummy));
3421: PetscCall(VecRestoreArrayRead(x, NULL));
3422: *a = NULL;
3423: PetscFunctionReturn(PETSC_SUCCESS);
3424: }
3426: /*@C
3427: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3428: processor's portion of the vector data. You MUST call `VecRestoreArray1dRead()`
3429: when you no longer need access to the array.
3431: Logically Collective
3433: Input Parameters:
3434: + x - the vector
3435: . m - first dimension of two dimensional array
3436: - mstart - first index you will use in first coordinate direction (often 0)
3438: Output Parameter:
3439: . a - location to put pointer to the array
3441: Level: developer
3443: Notes:
3444: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3445: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3446: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3448: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3450: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3451: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3452: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3453: @*/
3454: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3455: {
3456: PetscInt N;
3458: PetscFunctionBegin;
3460: PetscAssertPointer(a, 4);
3462: PetscCall(VecGetLocalSize(x, &N));
3463: PetscCheck(m == N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local array size %" PetscInt_FMT " does not match 1d array dimensions %" PetscInt_FMT, N, m);
3464: PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3465: *a -= mstart;
3466: PetscFunctionReturn(PETSC_SUCCESS);
3467: }
3469: /*@C
3470: VecRestoreArray1dRead - Restores a vector after `VecGetArray1dRead()` has been called.
3472: Logically Collective
3474: Input Parameters:
3475: + x - the vector
3476: . m - first dimension of two dimensional array
3477: . mstart - first index you will use in first coordinate direction (often 0)
3478: - a - location of pointer to array obtained from `VecGetArray1dRead()`
3480: Level: developer
3482: Notes:
3483: For regular PETSc vectors this routine does not involve any copies. For
3484: any special vectors that do not store local vector data in a contiguous
3485: array, this routine will copy the data back into the underlying
3486: vector data structure from the array obtained with `VecGetArray1dRead()`.
3488: This routine actually zeros out the `a` pointer.
3490: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3491: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3492: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3493: @*/
3494: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3495: {
3496: PetscFunctionBegin;
3499: PetscCall(VecRestoreArrayRead(x, NULL));
3500: *a = NULL;
3501: PetscFunctionReturn(PETSC_SUCCESS);
3502: }
3504: /*@C
3505: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3506: processor's portion of the vector data. You MUST call `VecRestoreArray3dRead()`
3507: when you no longer need access to the array.
3509: Logically Collective
3511: Input Parameters:
3512: + x - the vector
3513: . m - first dimension of three dimensional array
3514: . n - second dimension of three dimensional array
3515: . p - third dimension of three dimensional array
3516: . mstart - first index you will use in first coordinate direction (often 0)
3517: . nstart - first index in the second coordinate direction (often 0)
3518: - pstart - first index in the third coordinate direction (often 0)
3520: Output Parameter:
3521: . a - location to put pointer to the array
3523: Level: developer
3525: Notes:
3526: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3527: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3528: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3529: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3dRead()`.
3531: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3533: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3534: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3535: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3536: @*/
3537: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3538: {
3539: PetscInt i, N, j;
3540: const PetscScalar *aa;
3541: PetscScalar **b;
3543: PetscFunctionBegin;
3545: PetscAssertPointer(a, 8);
3547: PetscCall(VecGetLocalSize(x, &N));
3548: PetscCheck(m * n * p == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 3d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n, p);
3549: PetscCall(VecGetArrayRead(x, &aa));
3551: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3552: b = (PetscScalar **)((*a) + m);
3553: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3554: for (i = 0; i < m; i++)
3555: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset((PetscScalar *)aa, i * n * p + j * p - pstart);
3556: *a -= mstart;
3557: PetscFunctionReturn(PETSC_SUCCESS);
3558: }
3560: /*@C
3561: VecRestoreArray3dRead - Restores a vector after `VecGetArray3dRead()` has been called.
3563: Logically Collective
3565: Input Parameters:
3566: + x - the vector
3567: . m - first dimension of three dimensional array
3568: . n - second dimension of the three dimensional array
3569: . p - third dimension of the three dimensional array
3570: . mstart - first index you will use in first coordinate direction (often 0)
3571: . nstart - first index in the second coordinate direction (often 0)
3572: . pstart - first index in the third coordinate direction (often 0)
3573: - a - location of pointer to array obtained from `VecGetArray3dRead()`
3575: Level: developer
3577: Notes:
3578: For regular PETSc vectors this routine does not involve any copies. For
3579: any special vectors that do not store local vector data in a contiguous
3580: array, this routine will copy the data back into the underlying
3581: vector data structure from the array obtained with `VecGetArray()`.
3583: This routine actually zeros out the `a` pointer.
3585: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3586: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3587: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3588: @*/
3589: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3590: {
3591: void *dummy;
3593: PetscFunctionBegin;
3595: PetscAssertPointer(a, 8);
3597: dummy = (void *)(*a + mstart);
3598: PetscCall(PetscFree(dummy));
3599: PetscCall(VecRestoreArrayRead(x, NULL));
3600: *a = NULL;
3601: PetscFunctionReturn(PETSC_SUCCESS);
3602: }
3604: /*@C
3605: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3606: processor's portion of the vector data. You MUST call `VecRestoreArray4dRead()`
3607: when you no longer need access to the array.
3609: Logically Collective
3611: Input Parameters:
3612: + x - the vector
3613: . m - first dimension of four dimensional array
3614: . n - second dimension of four dimensional array
3615: . p - third dimension of four dimensional array
3616: . q - fourth dimension of four dimensional array
3617: . mstart - first index you will use in first coordinate direction (often 0)
3618: . nstart - first index in the second coordinate direction (often 0)
3619: . pstart - first index in the third coordinate direction (often 0)
3620: - qstart - first index in the fourth coordinate direction (often 0)
3622: Output Parameter:
3623: . a - location to put pointer to the array
3625: Level: beginner
3627: Notes:
3628: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3629: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3630: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3631: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3633: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3635: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3636: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3637: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3638: @*/
3639: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3640: {
3641: PetscInt i, N, j, k;
3642: const PetscScalar *aa;
3643: PetscScalar ***b, **c;
3645: PetscFunctionBegin;
3647: PetscAssertPointer(a, 10);
3649: PetscCall(VecGetLocalSize(x, &N));
3650: PetscCheck(m * n * p * q == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 4d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n, p, q);
3651: PetscCall(VecGetArrayRead(x, &aa));
3653: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3654: b = (PetscScalar ***)((*a) + m);
3655: c = (PetscScalar **)(b + m * n);
3656: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3657: for (i = 0; i < m; i++)
3658: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3659: for (i = 0; i < m; i++)
3660: for (j = 0; j < n; j++)
3661: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = (PetscScalar *)aa + i * n * p * q + j * p * q + k * q - qstart;
3662: *a -= mstart;
3663: PetscFunctionReturn(PETSC_SUCCESS);
3664: }
3666: /*@C
3667: VecRestoreArray4dRead - Restores a vector after `VecGetArray4d()` has been called.
3669: Logically Collective
3671: Input Parameters:
3672: + x - the vector
3673: . m - first dimension of four dimensional array
3674: . n - second dimension of the four dimensional array
3675: . p - third dimension of the four dimensional array
3676: . q - fourth dimension of the four dimensional array
3677: . mstart - first index you will use in first coordinate direction (often 0)
3678: . nstart - first index in the second coordinate direction (often 0)
3679: . pstart - first index in the third coordinate direction (often 0)
3680: . qstart - first index in the fourth coordinate direction (often 0)
3681: - a - location of pointer to array obtained from `VecGetArray4dRead()`
3683: Level: beginner
3685: Notes:
3686: For regular PETSc vectors this routine does not involve any copies. For
3687: any special vectors that do not store local vector data in a contiguous
3688: array, this routine will copy the data back into the underlying
3689: vector data structure from the array obtained with `VecGetArray()`.
3691: This routine actually zeros out the `a` pointer.
3693: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3694: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3695: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3696: @*/
3697: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3698: {
3699: void *dummy;
3701: PetscFunctionBegin;
3703: PetscAssertPointer(a, 10);
3705: dummy = (void *)(*a + mstart);
3706: PetscCall(PetscFree(dummy));
3707: PetscCall(VecRestoreArrayRead(x, NULL));
3708: *a = NULL;
3709: PetscFunctionReturn(PETSC_SUCCESS);
3710: }
3712: /*@
3713: VecLockGet - Get the current lock status of a vector
3715: Logically Collective
3717: Input Parameter:
3718: . x - the vector
3720: Output Parameter:
3721: . state - greater than zero indicates the vector is locked for read; less than zero indicates the vector is
3722: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3724: Level: advanced
3726: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3727: @*/
3728: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3729: {
3730: PetscFunctionBegin;
3732: PetscAssertPointer(state, 2);
3733: *state = x->lock;
3734: PetscFunctionReturn(PETSC_SUCCESS);
3735: }
3737: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3738: {
3739: PetscFunctionBegin;
3741: PetscAssertPointer(file, 2);
3742: PetscAssertPointer(func, 3);
3743: PetscAssertPointer(line, 4);
3744: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3745: {
3746: const int index = x->lockstack.currentsize - 1;
3748: *file = index < 0 ? NULL : x->lockstack.file[index];
3749: *func = index < 0 ? NULL : x->lockstack.function[index];
3750: *line = index < 0 ? 0 : x->lockstack.line[index];
3751: }
3752: #else
3753: *file = NULL;
3754: *func = NULL;
3755: *line = 0;
3756: #endif
3757: PetscFunctionReturn(PETSC_SUCCESS);
3758: }
3760: /*@
3761: VecLockReadPush - Push a read-only lock on a vector to prevent it from being written to
3763: Logically Collective
3765: Input Parameter:
3766: . x - the vector
3768: Level: intermediate
3770: Notes:
3771: If this is set then calls to `VecGetArray()` or `VecSetValues()` or any other routines that change the vectors values will generate an error.
3773: The call can be nested, i.e., called multiple times on the same vector, but each `VecLockReadPush()` has to have one matching
3774: `VecLockReadPop()`, which removes the latest read-only lock.
3776: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3777: @*/
3778: PetscErrorCode VecLockReadPush(Vec x)
3779: {
3780: PetscFunctionBegin;
3782: PetscCheck(x->lock++ >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector is already locked for exclusive write access but you want to read it");
3783: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3784: {
3785: const char *file, *func;
3786: int index, line;
3788: if ((index = petscstack.currentsize - 2) < 0) {
3789: // vec was locked "outside" of petsc, either in user-land or main. the error message will
3790: // now show this function as the culprit, but it will include the stacktrace
3791: file = "unknown user-file";
3792: func = "unknown_user_function";
3793: line = 0;
3794: } else {
3795: file = petscstack.file[index];
3796: func = petscstack.function[index];
3797: line = petscstack.line[index];
3798: }
3799: PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
3800: }
3801: #endif
3802: PetscFunctionReturn(PETSC_SUCCESS);
3803: }
3805: /*@
3806: VecLockReadPop - Pop a read-only lock from a vector
3808: Logically Collective
3810: Input Parameter:
3811: . x - the vector
3813: Level: intermediate
3815: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
3816: @*/
3817: PetscErrorCode VecLockReadPop(Vec x)
3818: {
3819: PetscFunctionBegin;
3821: PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
3822: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3823: {
3824: const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];
3826: PetscStackPop_Private(x->lockstack, previous);
3827: }
3828: #endif
3829: PetscFunctionReturn(PETSC_SUCCESS);
3830: }
3832: /*@
3833: VecLockWriteSet - Lock or unlock a vector for exclusive read/write access
3835: Logically Collective
3837: Input Parameters:
3838: + x - the vector
3839: - flg - `PETSC_TRUE` to lock the vector for exclusive read/write access; `PETSC_FALSE` to unlock it.
3841: Level: intermediate
3843: Notes:
3844: The function is useful in split-phase computations, which usually have a begin phase and an end phase.
3845: One can call `VecLockWriteSet`(x,`PETSC_TRUE`) in the begin phase to lock a vector for exclusive
3846: access, and call `VecLockWriteSet`(x,`PETSC_FALSE`) in the end phase to unlock the vector from exclusive
3847: access. In this way, one is ensured no other operations can access the vector in between. The code may like
3849: .vb
3850: VecGetArray(x,&xdata); // begin phase
3851: VecLockWriteSet(v,PETSC_TRUE);
3853: Other operations, which can not access x anymore (they can access xdata, of course)
3855: VecRestoreArray(x,&vdata); // end phase
3856: VecLockWriteSet(v,PETSC_FALSE);
3857: .ve
3859: The call can not be nested on the same vector, in other words, one can not call `VecLockWriteSet`(x,`PETSC_TRUE`)
3860: again before calling `VecLockWriteSet`(v,`PETSC_FALSE`).
3862: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
3863: @*/
3864: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
3865: {
3866: PetscFunctionBegin;
3868: if (flg) {
3869: PetscCheck(x->lock <= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector is already locked for read-only access but you want to write it");
3870: PetscCheck(x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector is already locked for exclusive write access but you want to write it");
3871: x->lock = -1;
3872: } else {
3873: PetscCheck(x->lock == -1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector is not locked for exclusive write access but you want to unlock it from that");
3874: x->lock = 0;
3875: }
3876: PetscFunctionReturn(PETSC_SUCCESS);
3877: }