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: .vb
94: val = (x,y) = y^H x,
95: .ve
96: where y^H denotes the conjugate transpose of y. Note that this corresponds to the usual "mathematicians" complex
97: inner product where the SECOND argument gets the complex conjugate. Since the `BLASdot()` complex conjugates the first
98: first argument we call the `BLASdot()` with the arguments reversed.
100: Use `VecTDot()` for the indefinite form
101: .vb
102: val = (x,y) = y^T x,
103: .ve
104: where y^T denotes the transpose of y.
106: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDotRealPart()`
107: @*/
108: PetscErrorCode VecDot(Vec x, Vec y, PetscScalar *val)
109: {
110: PetscFunctionBegin;
113: PetscAssertPointer(val, 3);
116: PetscCheckSameTypeAndComm(x, 1, y, 2);
117: VecCheckSameSize(x, 1, y, 2);
118: VecCheckAssembled(x);
119: VecCheckAssembled(y);
121: PetscCall(VecLockReadPush(x));
122: PetscCall(VecLockReadPush(y));
123: PetscCall(PetscLogEventBegin(VEC_Dot, x, y, 0, 0));
124: PetscUseTypeMethod(x, dot, y, val);
125: PetscCall(PetscLogEventEnd(VEC_Dot, x, y, 0, 0));
126: PetscCall(VecLockReadPop(x));
127: PetscCall(VecLockReadPop(y));
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: /*@
132: VecDotRealPart - Computes the real part of the vector dot product.
134: Collective
136: Input Parameters:
137: + x - first vector
138: - y - second vector
140: Output Parameter:
141: . val - the real part of the dot product;
143: Level: intermediate
145: Notes for Users of Complex Numbers:
146: See `VecDot()` for more details on the definition of the dot product for complex numbers
148: For real numbers this returns the same value as `VecDot()`
150: 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
151: the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)
153: Developer Notes:
154: This is not currently optimized to compute only the real part of the dot product.
156: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDot()`, `VecDotNorm2()`
157: @*/
158: PetscErrorCode VecDotRealPart(Vec x, Vec y, PetscReal *val)
159: {
160: PetscScalar fdot;
162: PetscFunctionBegin;
163: PetscCall(VecDot(x, y, &fdot));
164: *val = PetscRealPart(fdot);
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: /*@
169: VecNorm - Computes the vector norm.
171: Collective
173: Input Parameters:
174: + x - the vector
175: - type - the type of the norm requested
177: Output Parameter:
178: . val - the norm
180: Level: intermediate
182: Notes:
183: See `NormType` for descriptions of each norm.
185: For complex numbers `NORM_1` will return the traditional 1 norm of the 2 norm of the complex
186: numbers; that is the 1 norm of the absolute values of the complex entries. In PETSc 3.6 and
187: earlier releases it returned the 1 norm of the 1 norm of the complex entries (what is
188: returned by the BLAS routine `asum()`). Both are valid norms but most people expect the former.
190: This routine stashes the computed norm value, repeated calls before the vector entries are
191: changed are then rapid since the precomputed value is immediately available. Certain vector
192: operations such as `VecSet()` store the norms so the value is immediately available and does
193: not need to be explicitly computed. `VecScale()` updates any stashed norm values, thus calls
194: after `VecScale()` do not need to explicitly recompute the norm.
196: .seealso: [](ch_vectors), `Vec`, `NormType`, `VecDot()`, `VecTDot()`, `VecDotBegin()`, `VecDotEnd()`, `VecNormAvailable()`,
197: `VecNormBegin()`, `VecNormEnd()`, `NormType()`
198: @*/
199: PetscErrorCode VecNorm(Vec x, NormType type, PetscReal *val)
200: {
201: PetscBool flg = PETSC_TRUE;
203: PetscFunctionBegin;
206: VecCheckAssembled(x);
208: PetscAssertPointer(val, 3);
210: PetscCall(VecNormAvailable(x, type, &flg, val));
211: // check that all MPI processes call this routine together and have same availability
212: if (PetscDefined(USE_DEBUG)) {
213: PetscMPIInt b0 = (PetscMPIInt)flg, b1[2], b2[2];
214: b1[0] = -b0;
215: b1[1] = b0;
216: PetscCallMPI(MPIU_Allreduce(b1, b2, 2, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
217: 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]);
218: if (flg) {
219: PetscReal b1[2], b2[2];
220: b1[0] = -(*val);
221: b1[1] = *val;
222: PetscCallMPI(MPIU_Allreduce(b1, b2, 2, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)x)));
223: 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);
224: }
225: }
226: if (flg) PetscFunctionReturn(PETSC_SUCCESS);
228: PetscCall(VecLockReadPush(x));
229: PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
230: PetscUseTypeMethod(x, norm, type, val);
231: PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
232: PetscCall(VecLockReadPop(x));
234: if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: /*@
239: VecNormAvailable - Returns the vector norm if it is already known. That is, it has been previously computed and cached in the vector
241: Not Collective
243: Input Parameters:
244: + x - the vector
245: - 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
246: `NORM_1_AND_2`, which computes both norms and stores them
247: in a two element array.
249: Output Parameters:
250: + available - `PETSC_TRUE` if the val returned is valid
251: - val - the norm
253: Level: intermediate
255: Developer Notes:
256: `PETSC_HAVE_SLOW_BLAS_NORM2` will cause a C (loop unrolled) version of the norm to be used, rather
257: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
258: (as opposed to vendor provided) because the FORTRAN BLAS `NRM2()` routine is very slow.
260: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
261: `VecNormBegin()`, `VecNormEnd()`
262: @*/
263: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
264: {
265: PetscFunctionBegin;
268: PetscAssertPointer(available, 3);
269: PetscAssertPointer(val, 4);
271: if (type == NORM_1_AND_2) {
272: *available = PETSC_FALSE;
273: } else {
274: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
275: }
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /*@
280: VecNormalize - Normalizes a vector by its 2-norm.
282: Collective
284: Input Parameter:
285: . x - the vector
287: Output Parameter:
288: . val - the vector norm before normalization. May be `NULL` if the value is not needed.
290: Level: intermediate
292: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
293: @*/
294: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
295: {
296: PetscReal norm;
298: PetscFunctionBegin;
301: PetscCall(VecSetErrorIfLocked(x, 1));
302: if (val) PetscAssertPointer(val, 2);
303: PetscCall(PetscLogEventBegin(VEC_Normalize, x, 0, 0, 0));
304: PetscCall(VecNorm(x, NORM_2, &norm));
305: if (norm == 0.0) PetscCall(PetscInfo(x, "Vector of zero norm can not be normalized; Returning only the zero norm\n"));
306: else if (PetscIsInfOrNanReal(norm)) PetscCall(PetscInfo(x, "Vector with Inf or Nan norm can not be normalized; Returning only the norm\n"));
307: else {
308: PetscScalar s = 1.0 / norm;
309: PetscCall(VecScale(x, s));
310: }
311: PetscCall(PetscLogEventEnd(VEC_Normalize, x, 0, 0, 0));
312: if (val) *val = norm;
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: /*@
317: VecMax - Determines the vector component with maximum real part and its location.
319: Collective
321: Input Parameter:
322: . x - the vector
324: Output Parameters:
325: + p - the index of `val` (pass `NULL` if you don't want this) in the vector
326: - val - the maximum component
328: Level: intermediate
330: Notes:
331: Returns the value `PETSC_MIN_REAL` and negative `p` if the vector is of length 0.
333: Returns the smallest index with the maximum value
335: Developer Note:
336: The Nag Fortran compiler does not like the symbol name VecMax
338: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
339: @*/
340: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
341: {
342: PetscFunctionBegin;
345: VecCheckAssembled(x);
346: if (p) PetscAssertPointer(p, 2);
347: PetscAssertPointer(val, 3);
348: PetscCall(VecLockReadPush(x));
349: PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
350: PetscUseTypeMethod(x, max, p, val);
351: PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
352: PetscCall(VecLockReadPop(x));
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: /*@
357: VecMin - Determines the vector component with minimum real part and its location.
359: Collective
361: Input Parameter:
362: . x - the vector
364: Output Parameters:
365: + p - the index of `val` (pass `NULL` if you don't want this location) in the vector
366: - val - the minimum component
368: Level: intermediate
370: Notes:
371: Returns the value `PETSC_MAX_REAL` and negative `p` if the vector is of length 0.
373: This returns the smallest index with the minimum value
375: Developer Note:
376: The Nag Fortran compiler does not like the symbol name VecMin
378: .seealso: [](ch_vectors), `Vec`, `VecMax()`
379: @*/
380: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
381: {
382: PetscFunctionBegin;
385: VecCheckAssembled(x);
386: if (p) PetscAssertPointer(p, 2);
387: PetscAssertPointer(val, 3);
388: PetscCall(VecLockReadPush(x));
389: PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
390: PetscUseTypeMethod(x, min, p, val);
391: PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
392: PetscCall(VecLockReadPop(x));
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: /*@
397: VecTDot - Computes an indefinite vector dot product. That is, this
398: routine does NOT use the complex conjugate.
400: Collective
402: Input Parameters:
403: + x - first vector
404: - y - second vector
406: Output Parameter:
407: . val - the dot product
409: Level: intermediate
411: Notes for Users of Complex Numbers:
412: For complex vectors, `VecTDot()` computes the indefinite form
413: .vb
414: val = (x,y) = y^T x,
415: .ve
416: where y^T denotes the transpose of y.
418: Use `VecDot()` for the inner product
419: .vb
420: val = (x,y) = y^H x,
421: .ve
422: where y^H denotes the conjugate transpose of y.
424: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
425: @*/
426: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
427: {
428: PetscFunctionBegin;
431: PetscAssertPointer(val, 3);
434: PetscCheckSameTypeAndComm(x, 1, y, 2);
435: VecCheckSameSize(x, 1, y, 2);
436: VecCheckAssembled(x);
437: VecCheckAssembled(y);
439: PetscCall(VecLockReadPush(x));
440: PetscCall(VecLockReadPush(y));
441: PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
442: PetscUseTypeMethod(x, tdot, y, val);
443: PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
444: PetscCall(VecLockReadPop(x));
445: PetscCall(VecLockReadPop(y));
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: PetscErrorCode VecScaleAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
450: {
451: PetscReal norms[4];
452: PetscBool flgs[4];
453: PetscScalar one = 1.0;
455: PetscFunctionBegin;
458: VecCheckAssembled(x);
459: PetscCall(VecSetErrorIfLocked(x, 1));
461: if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);
463: /* get current stashed norms */
464: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
466: PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
467: VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
468: PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));
470: PetscCall(PetscObjectStateIncrease((PetscObject)x));
471: /* put the scaled stashed norms back into the Vec */
472: for (PetscInt i = 0; i < 4; i++) {
473: PetscReal ar = PetscAbsScalar(alpha);
474: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
475: }
476: PetscFunctionReturn(PETSC_SUCCESS);
477: }
479: /*@
480: VecScale - Scales a vector.
482: Logically Collective
484: Input Parameters:
485: + x - the vector
486: - alpha - the scalar
488: Level: intermediate
490: Note:
491: For a vector with n components, `VecScale()` computes x[i] = alpha * x[i], for i=1,...,n.
493: .seealso: [](ch_vectors), `Vec`, `VecSet()`
494: @*/
495: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
496: {
497: PetscFunctionBegin;
498: PetscCall(VecScaleAsync_Private(x, alpha, NULL));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
503: {
504: PetscFunctionBegin;
507: VecCheckAssembled(x);
509: PetscCall(VecSetErrorIfLocked(x, 1));
511: if (alpha == 0) {
512: PetscReal norm;
513: PetscBool set;
515: PetscCall(VecNormAvailable(x, NORM_2, &set, &norm));
516: if (set == PETSC_TRUE && norm == 0) PetscFunctionReturn(PETSC_SUCCESS);
517: }
518: PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
519: VecMethodDispatch(x, dctx, VecAsyncFnName(Set), set, (Vec, PetscScalar, PetscDeviceContext), alpha);
520: PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
521: PetscCall(PetscObjectStateIncrease((PetscObject)x));
523: /* norms can be simply set (if |alpha|*N not too large) */
524: {
525: PetscReal val = PetscAbsScalar(alpha);
526: const PetscInt N = x->map->N;
528: if (N == 0) {
529: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0));
530: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
531: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
532: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
533: } else if (val > PETSC_MAX_REAL / N) {
534: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
535: } else {
536: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
537: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
538: val *= PetscSqrtReal((PetscReal)N);
539: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
540: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
541: }
542: }
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: /*@
547: VecSet - Sets all components of a vector to a single scalar value.
549: Logically Collective
551: Input Parameters:
552: + x - the vector
553: - alpha - the scalar
555: Level: beginner
557: Notes:
558: For a vector of dimension n, `VecSet()` sets x[i] = alpha, for i=1,...,n,
559: so that all vector entries then equal the identical
560: scalar value, `alpha`. Use the more general routine
561: `VecSetValues()` to set different vector entries.
563: You CANNOT call this after you have called `VecSetValues()` but before you call
564: `VecAssemblyBegin()`
566: If `alpha` is zero and the norm of the vector is known to be zero then this skips the unneeded zeroing process
568: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
569: @*/
570: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
571: {
572: PetscFunctionBegin;
573: PetscCall(VecSetAsync_Private(x, alpha, NULL));
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: PetscErrorCode VecAXPYAsync_Private(Vec y, PetscScalar alpha, Vec x, PetscDeviceContext dctx)
578: {
579: PetscFunctionBegin;
584: PetscCheckSameTypeAndComm(x, 3, y, 1);
585: VecCheckSameSize(x, 3, y, 1);
586: VecCheckAssembled(x);
587: VecCheckAssembled(y);
589: if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
590: PetscCall(VecSetErrorIfLocked(y, 1));
591: if (x == y) {
592: PetscCall(VecScale(y, alpha + 1.0));
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
595: PetscCall(VecLockReadPush(x));
596: PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
597: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPY), axpy, (Vec, PetscScalar, Vec, PetscDeviceContext), alpha, x);
598: PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
599: PetscCall(VecLockReadPop(x));
600: PetscCall(PetscObjectStateIncrease((PetscObject)y));
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }
603: /*@
604: VecAXPY - Computes `y = alpha x + y`.
606: Logically Collective
608: Input Parameters:
609: + alpha - the scalar
610: . x - vector scale by `alpha`
611: - y - vector accumulated into
613: Output Parameter:
614: . y - output vector
616: Level: intermediate
618: Notes:
619: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
620: .vb
621: VecAXPY(y,alpha,x) y = alpha x + y
622: VecAYPX(y,beta,x) y = x + beta y
623: VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
624: VecWAXPY(w,alpha,x,y) w = alpha x + y
625: VecAXPBYPCZ(z,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
626: VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
627: .ve
629: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
630: @*/
631: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
632: {
633: PetscFunctionBegin;
634: PetscCall(VecAXPYAsync_Private(y, alpha, x, NULL));
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: PetscErrorCode VecAYPXAsync_Private(Vec y, PetscScalar beta, Vec x, PetscDeviceContext dctx)
639: {
640: PetscFunctionBegin;
645: PetscCheckSameTypeAndComm(x, 3, y, 1);
646: VecCheckSameSize(x, 1, y, 3);
647: VecCheckAssembled(x);
648: VecCheckAssembled(y);
650: PetscCall(VecSetErrorIfLocked(y, 1));
651: if (x == y) {
652: PetscCall(VecScale(y, beta + 1.0));
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
655: PetscCall(VecLockReadPush(x));
656: if (beta == (PetscScalar)0.0) {
657: PetscCall(VecCopy(x, y));
658: } else {
659: PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
660: VecMethodDispatch(y, dctx, VecAsyncFnName(AYPX), aypx, (Vec, PetscScalar, Vec, PetscDeviceContext), beta, x);
661: PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
662: PetscCall(PetscObjectStateIncrease((PetscObject)y));
663: }
664: PetscCall(VecLockReadPop(x));
665: PetscFunctionReturn(PETSC_SUCCESS);
666: }
668: /*@
669: VecAYPX - Computes `y = x + beta y`.
671: Logically Collective
673: Input Parameters:
674: + beta - the scalar
675: . x - the unscaled vector
676: - y - the vector to be scaled
678: Output Parameter:
679: . y - output vector
681: Level: intermediate
683: Developer Notes:
684: The implementation is optimized for `beta` of -1.0, 0.0, and 1.0
686: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
687: @*/
688: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
689: {
690: PetscFunctionBegin;
691: PetscCall(VecAYPXAsync_Private(y, beta, x, NULL));
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: PetscErrorCode VecAXPBYAsync_Private(Vec y, PetscScalar alpha, PetscScalar beta, Vec x, PetscDeviceContext dctx)
696: {
697: PetscFunctionBegin;
702: PetscCheckSameTypeAndComm(x, 4, y, 1);
703: VecCheckSameSize(y, 1, x, 4);
704: VecCheckAssembled(x);
705: VecCheckAssembled(y);
708: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
709: if (x == y) {
710: PetscCall(VecScale(y, alpha + beta));
711: PetscFunctionReturn(PETSC_SUCCESS);
712: }
714: PetscCall(VecSetErrorIfLocked(y, 1));
715: PetscCall(VecLockReadPush(x));
716: PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
717: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPBY), axpby, (Vec, PetscScalar, PetscScalar, Vec, PetscDeviceContext), alpha, beta, x);
718: PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
719: PetscCall(PetscObjectStateIncrease((PetscObject)y));
720: PetscCall(VecLockReadPop(x));
721: PetscFunctionReturn(PETSC_SUCCESS);
722: }
724: /*@
725: VecAXPBY - Computes `y = alpha x + beta y`.
727: Logically Collective
729: Input Parameters:
730: + alpha - first scalar
731: . beta - second scalar
732: . x - the first scaled vector
733: - y - the second scaled vector
735: Output Parameter:
736: . y - output vector
738: Level: intermediate
740: Developer Notes:
741: The implementation is optimized for `alpha` and/or `beta` values of 0.0 and 1.0
743: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
744: @*/
745: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
746: {
747: PetscFunctionBegin;
748: PetscCall(VecAXPBYAsync_Private(y, alpha, beta, x, NULL));
749: PetscFunctionReturn(PETSC_SUCCESS);
750: }
752: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
753: {
754: PetscFunctionBegin;
761: PetscCheckSameTypeAndComm(x, 5, y, 6);
762: PetscCheckSameTypeAndComm(x, 5, z, 1);
763: VecCheckSameSize(x, 5, y, 6);
764: VecCheckSameSize(x, 5, z, 1);
765: PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
766: PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
767: VecCheckAssembled(x);
768: VecCheckAssembled(y);
769: VecCheckAssembled(z);
773: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
775: PetscCall(VecSetErrorIfLocked(z, 1));
776: PetscCall(VecLockReadPush(x));
777: PetscCall(VecLockReadPush(y));
778: PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
779: VecMethodDispatch(z, dctx, VecAsyncFnName(AXPBYPCZ), axpbypcz, (Vec, PetscScalar, PetscScalar, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, beta, gamma, x, y);
780: PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
781: PetscCall(PetscObjectStateIncrease((PetscObject)z));
782: PetscCall(VecLockReadPop(x));
783: PetscCall(VecLockReadPop(y));
784: PetscFunctionReturn(PETSC_SUCCESS);
785: }
786: /*@
787: VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`
789: Logically Collective
791: Input Parameters:
792: + alpha - first scalar
793: . beta - second scalar
794: . gamma - third scalar
795: . x - first vector
796: . y - second vector
797: - z - third vector
799: Output Parameter:
800: . z - output vector
802: Level: intermediate
804: Note:
805: `x`, `y` and `z` must be different vectors
807: Developer Notes:
808: The implementation is optimized for `alpha` of 1.0 and `gamma` of 1.0 or 0.0
810: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
811: @*/
812: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
813: {
814: PetscFunctionBegin;
815: PetscCall(VecAXPBYPCZAsync_Private(z, alpha, beta, gamma, x, y, NULL));
816: PetscFunctionReturn(PETSC_SUCCESS);
817: }
819: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
820: {
821: PetscFunctionBegin;
828: PetscCheckSameTypeAndComm(x, 3, y, 4);
829: PetscCheckSameTypeAndComm(y, 4, w, 1);
830: VecCheckSameSize(x, 3, y, 4);
831: VecCheckSameSize(x, 3, w, 1);
832: PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
833: PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
834: VecCheckAssembled(x);
835: VecCheckAssembled(y);
837: PetscCall(VecSetErrorIfLocked(w, 1));
839: PetscCall(VecLockReadPush(x));
840: PetscCall(VecLockReadPush(y));
841: if (alpha == (PetscScalar)0.0) {
842: PetscCall(VecCopyAsync_Private(y, w, dctx));
843: } else {
844: PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
845: VecMethodDispatch(w, dctx, VecAsyncFnName(WAXPY), waxpy, (Vec, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, x, y);
846: PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
847: PetscCall(PetscObjectStateIncrease((PetscObject)w));
848: }
849: PetscCall(VecLockReadPop(x));
850: PetscCall(VecLockReadPop(y));
851: PetscFunctionReturn(PETSC_SUCCESS);
852: }
854: /*@
855: VecWAXPY - Computes `w = alpha x + y`.
857: Logically Collective
859: Input Parameters:
860: + alpha - the scalar
861: . x - first vector, multiplied by `alpha`
862: - y - second vector
864: Output Parameter:
865: . w - the result
867: Level: intermediate
869: Note:
870: `w` cannot be either `x` or `y`, but `x` and `y` can be the same
872: Developer Notes:
873: The implementation is optimized for alpha of -1.0, 0.0, and 1.0
875: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
876: @*/
877: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
878: {
879: PetscFunctionBegin;
880: PetscCall(VecWAXPYAsync_Private(w, alpha, x, y, NULL));
881: PetscFunctionReturn(PETSC_SUCCESS);
882: }
884: /*@
885: VecSetValues - Inserts or adds values into certain locations of a vector.
887: Not Collective
889: Input Parameters:
890: + x - vector to insert in
891: . ni - number of elements to add
892: . ix - indices where to add
893: . y - array of values
894: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries
896: Level: beginner
898: Notes:
899: .vb
900: `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
901: .ve
903: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
904: options cannot be mixed without intervening calls to the assembly
905: routines.
907: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
908: MUST be called after all calls to `VecSetValues()` have been completed.
910: VecSetValues() uses 0-based indices in Fortran as well as in C.
912: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
913: negative indices may be passed in ix. These rows are
914: simply ignored. This allows easily inserting element load matrices
915: with homogeneous Dirichlet boundary conditions that you don't want represented
916: in the vector.
918: Fortran Note:
919: If any of `ix` and `y` are scalars pass them using, for example,
920: .vb
921: call VecSetValues(mat, one, [ix], [y], INSERT_VALUES, ierr)
922: .ve
924: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
925: `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
926: @*/
927: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
928: {
929: PetscFunctionBeginHot;
931: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
932: PetscAssertPointer(ix, 3);
933: PetscAssertPointer(y, 4);
936: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
937: PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
938: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
939: PetscCall(PetscObjectStateIncrease((PetscObject)x));
940: PetscFunctionReturn(PETSC_SUCCESS);
941: }
943: /*@
944: VecGetValues - Gets values from certain locations of a vector. Currently
945: can only get values on the same processor on which they are owned
947: Not Collective
949: Input Parameters:
950: + x - vector to get values from
951: . ni - number of elements to get
952: - ix - indices where to get them from (in global 1d numbering)
954: Output Parameter:
955: . y - array of values, must be passed in with a length of `ni`
957: Level: beginner
959: Notes:
960: The user provides the allocated array y; it is NOT allocated in this routine
962: `VecGetValues()` gets y[i] = x[ix[i]], for i=0,...,ni-1.
964: `VecAssemblyBegin()` and `VecAssemblyEnd()` MUST be called before calling this if `VecSetValues()` or related routine has been called
966: VecGetValues() uses 0-based indices in Fortran as well as in C.
968: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
969: negative indices may be passed in ix. These rows are
970: simply ignored.
972: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
973: @*/
974: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
975: {
976: PetscFunctionBegin;
978: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
979: PetscAssertPointer(ix, 3);
980: PetscAssertPointer(y, 4);
982: VecCheckAssembled(x);
983: PetscUseTypeMethod(x, getvalues, ni, ix, y);
984: PetscFunctionReturn(PETSC_SUCCESS);
985: }
987: /*@
988: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
990: Not Collective
992: Input Parameters:
993: + x - vector to insert in
994: . ni - number of blocks to add
995: . ix - indices where to add in block count, rather than element count
996: . y - array of values
997: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries
999: Level: intermediate
1001: Notes:
1002: `VecSetValuesBlocked()` sets x[bs*ix[i]+j] = y[bs*i+j],
1003: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
1005: Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
1006: options cannot be mixed without intervening calls to the assembly
1007: routines.
1009: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1010: MUST be called after all calls to `VecSetValuesBlocked()` have been completed.
1012: `VecSetValuesBlocked()` uses 0-based indices in Fortran as well as in C.
1014: Negative indices may be passed in ix, these rows are
1015: simply ignored. This allows easily inserting element load matrices
1016: with homogeneous Dirichlet boundary conditions that you don't want represented
1017: in the vector.
1019: Fortran Note:
1020: If any of `ix` and `y` are scalars pass them using, for example,
1021: .vb
1022: call VecSetValuesBlocked(mat, one, [ix], [y], INSERT_VALUES, ierr)
1023: .ve
1025: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
1026: `VecSetValues()`
1027: @*/
1028: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1029: {
1030: PetscFunctionBeginHot;
1032: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1033: PetscAssertPointer(ix, 3);
1034: PetscAssertPointer(y, 4);
1037: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1038: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1039: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1040: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1041: PetscFunctionReturn(PETSC_SUCCESS);
1042: }
1044: /*@
1045: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1046: using a local ordering of the nodes.
1048: Not Collective
1050: Input Parameters:
1051: + x - vector to insert in
1052: . ni - number of elements to add
1053: . ix - indices where to add
1054: . y - array of values
1055: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1057: Level: intermediate
1059: Notes:
1060: `VecSetValuesLocal()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
1062: Calls to `VecSetValuesLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1063: options cannot be mixed without intervening calls to the assembly
1064: routines.
1066: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1067: MUST be called after all calls to `VecSetValuesLocal()` have been completed.
1069: `VecSetValuesLocal()` uses 0-based indices in Fortran as well as in C.
1071: Fortran Note:
1072: If any of `ix` and `y` are scalars pass them using, for example,
1073: .vb
1074: call VecSetValuesLocal(mat, one, [ix], [y], INSERT_VALUES, ierr)
1075: .ve
1077: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1078: `VecSetValuesBlockedLocal()`
1079: @*/
1080: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1081: {
1082: PetscInt lixp[128], *lix = lixp;
1084: PetscFunctionBeginHot;
1086: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1087: PetscAssertPointer(ix, 3);
1088: PetscAssertPointer(y, 4);
1091: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1092: if (!x->ops->setvalueslocal) {
1093: if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
1094: if (x->map->mapping) {
1095: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1096: PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1097: PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1098: if (ni > 128) PetscCall(PetscFree(lix));
1099: } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1100: } else PetscUseTypeMethod(x, setvalueslocal, ni, ix, y, iora);
1101: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1102: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1103: PetscFunctionReturn(PETSC_SUCCESS);
1104: }
1106: /*@
1107: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1108: using a local ordering of the nodes.
1110: Not Collective
1112: Input Parameters:
1113: + x - vector to insert in
1114: . ni - number of blocks to add
1115: . ix - indices where to add in block count, not element count
1116: . y - array of values
1117: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1119: Level: intermediate
1121: Notes:
1122: `VecSetValuesBlockedLocal()` sets x[bs*ix[i]+j] = y[bs*i+j],
1123: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with `VecSetBlockSize()`.
1125: Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1126: options cannot be mixed without intervening calls to the assembly
1127: routines.
1129: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1130: MUST be called after all calls to `VecSetValuesBlockedLocal()` have been completed.
1132: `VecSetValuesBlockedLocal()` uses 0-based indices in Fortran as well as in C.
1134: Fortran Note:
1135: If any of `ix` and `y` are scalars pass them using, for example,
1136: .vb
1137: call VecSetValuesBlockedLocal(mat, one, [ix], [y], INSERT_VALUES, ierr)
1138: .ve
1140: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1141: `VecSetLocalToGlobalMapping()`
1142: @*/
1143: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1144: {
1145: PetscInt lixp[128], *lix = lixp;
1147: PetscFunctionBeginHot;
1149: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1150: PetscAssertPointer(ix, 3);
1151: PetscAssertPointer(y, 4);
1153: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1154: if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
1155: if (x->map->mapping) {
1156: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1157: PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1158: PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1159: if (ni > 128) PetscCall(PetscFree(lix));
1160: } else {
1161: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1162: }
1163: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1164: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1165: PetscFunctionReturn(PETSC_SUCCESS);
1166: }
1168: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1169: {
1170: PetscFunctionBegin;
1173: VecCheckAssembled(x);
1175: if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1176: PetscAssertPointer(y, 3);
1177: for (PetscInt i = 0; i < nv; ++i) {
1180: PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1181: VecCheckSameSize(x, 1, y[i], 3);
1182: VecCheckAssembled(y[i]);
1183: PetscCall(VecLockReadPush(y[i]));
1184: }
1185: PetscAssertPointer(result, 4);
1188: PetscCall(VecLockReadPush(x));
1189: PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1190: PetscCall((*mxdot)(x, nv, y, result));
1191: PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1192: PetscCall(VecLockReadPop(x));
1193: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1194: PetscFunctionReturn(PETSC_SUCCESS);
1195: }
1197: /*@
1198: VecMTDot - Computes indefinite vector multiple dot products.
1199: That is, it does NOT use the complex conjugate.
1201: Collective
1203: Input Parameters:
1204: + x - one vector
1205: . nv - number of vectors
1206: - y - array of vectors. Note that vectors are pointers
1208: Output Parameter:
1209: . val - array of the dot products
1211: Level: intermediate
1213: Notes for Users of Complex Numbers:
1214: For complex vectors, `VecMTDot()` computes the indefinite form
1215: .vb
1216: val = (x,y) = y^T x,
1217: .ve
1218: where y^T denotes the transpose of y.
1220: Use `VecMDot()` for the inner product
1221: .vb
1222: val = (x,y) = y^H x,
1223: .ve
1224: where y^H denotes the conjugate transpose of y.
1226: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1227: @*/
1228: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1229: {
1230: PetscFunctionBegin;
1232: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1233: PetscFunctionReturn(PETSC_SUCCESS);
1234: }
1236: /*@
1237: VecMDot - Computes multiple vector dot products.
1239: Collective
1241: Input Parameters:
1242: + x - one vector
1243: . nv - number of vectors
1244: - y - array of vectors.
1246: Output Parameter:
1247: . val - array of the dot products (does not allocate the array)
1249: Level: intermediate
1251: Notes for Users of Complex Numbers:
1252: For complex vectors, `VecMDot()` computes
1253: .vb
1254: val = (x,y) = y^H x,
1255: .ve
1256: where y^H denotes the conjugate transpose of y.
1258: Use `VecMTDot()` for the indefinite form
1259: .vb
1260: val = (x,y) = y^T x,
1261: .ve
1262: where y^T denotes the transpose of y.
1264: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`
1265: @*/
1266: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1267: {
1268: PetscFunctionBegin;
1270: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1271: PetscFunctionReturn(PETSC_SUCCESS);
1272: }
1274: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
1275: {
1276: PetscFunctionBegin;
1278: VecCheckAssembled(y);
1280: PetscCall(VecSetErrorIfLocked(y, 1));
1281: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1282: if (nv) {
1283: PetscInt zeros = 0;
1285: PetscAssertPointer(alpha, 3);
1286: PetscAssertPointer(x, 4);
1287: for (PetscInt i = 0; i < nv; ++i) {
1291: PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1292: VecCheckSameSize(y, 1, x[i], 4);
1293: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1294: VecCheckAssembled(x[i]);
1295: PetscCall(VecLockReadPush(x[i]));
1296: zeros += alpha[i] == (PetscScalar)0.0;
1297: }
1299: if (zeros < nv) {
1300: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1301: VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1302: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1303: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1304: }
1306: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1307: }
1308: PetscFunctionReturn(PETSC_SUCCESS);
1309: }
1311: /*@
1312: VecMAXPY - Computes `y = y + sum alpha[i] x[i]`
1314: Logically Collective
1316: Input Parameters:
1317: + nv - number of scalars and x-vectors
1318: . alpha - array of scalars
1319: . y - one vector
1320: - x - array of vectors
1322: Level: intermediate
1324: Note:
1325: `y` cannot be any of the `x` vectors
1327: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`,`VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1328: @*/
1329: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1330: {
1331: PetscFunctionBegin;
1332: PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1333: PetscFunctionReturn(PETSC_SUCCESS);
1334: }
1336: /*@
1337: VecMAXPBY - Computes `y = beta y + sum alpha[i] x[i]`
1339: Logically Collective
1341: Input Parameters:
1342: + nv - number of scalars and x-vectors
1343: . alpha - array of scalars
1344: . beta - scalar
1345: . y - one vector
1346: - x - array of vectors
1348: Level: intermediate
1350: Note:
1351: `y` cannot be any of the `x` vectors.
1353: Developer Notes:
1354: This is a convenience routine, but implementations might be able to optimize it, for example, when `beta` is zero.
1356: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1357: @*/
1358: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
1359: {
1360: PetscFunctionBegin;
1362: VecCheckAssembled(y);
1364: PetscCall(VecSetErrorIfLocked(y, 1));
1365: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1368: if (y->ops->maxpby) {
1369: PetscInt zeros = 0;
1371: if (nv) {
1372: PetscAssertPointer(alpha, 3);
1373: PetscAssertPointer(x, 5);
1374: }
1376: for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
1380: PetscCheckSameTypeAndComm(y, 1, x[i], 5);
1381: VecCheckSameSize(y, 1, x[i], 5);
1382: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1383: VecCheckAssembled(x[i]);
1384: PetscCall(VecLockReadPush(x[i]));
1385: zeros += alpha[i] == (PetscScalar)0.0;
1386: }
1388: if (zeros < nv) { // has nonzero alpha
1389: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1390: PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1391: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1392: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1393: } else {
1394: PetscCall(VecScale(y, beta));
1395: }
1397: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1398: } else { // no maxpby
1399: if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1400: else PetscCall(VecScale(y, beta));
1401: PetscCall(VecMAXPY(y, nv, alpha, x));
1402: }
1403: PetscFunctionReturn(PETSC_SUCCESS);
1404: }
1406: /*@
1407: VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
1408: in the order they appear in the array. The concatenated vector resides on the same
1409: communicator and is the same type as the source vectors.
1411: Collective
1413: Input Parameters:
1414: + nx - number of vectors to be concatenated
1415: - X - array containing the vectors to be concatenated in the order of concatenation
1417: Output Parameters:
1418: + Y - concatenated vector
1419: - x_is - array of index sets corresponding to the concatenated components of `Y` (pass `NULL` if not needed)
1421: Level: advanced
1423: Notes:
1424: Concatenation is similar to the functionality of a `VECNEST` object; they both represent combination of
1425: different vector spaces. However, concatenated vectors do not store any information about their
1426: sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1427: manipulation of data in the concatenated vector that corresponds to the original components at creation.
1429: This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1430: has to operate on combined vector spaces and cannot utilize `VECNEST` objects due to incompatibility with
1431: bound projections.
1433: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1434: @*/
1435: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1436: {
1437: MPI_Comm comm;
1438: VecType vec_type;
1439: Vec Ytmp, Xtmp;
1440: IS *is_tmp;
1441: PetscInt i, shift = 0, Xnl, Xng, Xbegin;
1443: PetscFunctionBegin;
1447: PetscAssertPointer(Y, 3);
1449: if ((*X)->ops->concatenate) {
1450: /* use the dedicated concatenation function if available */
1451: PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1452: } else {
1453: /* loop over vectors and start creating IS */
1454: comm = PetscObjectComm((PetscObject)*X);
1455: PetscCall(VecGetType(*X, &vec_type));
1456: PetscCall(PetscMalloc1(nx, &is_tmp));
1457: for (i = 0; i < nx; i++) {
1458: PetscCall(VecGetSize(X[i], &Xng));
1459: PetscCall(VecGetLocalSize(X[i], &Xnl));
1460: PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1461: PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1462: shift += Xng;
1463: }
1464: /* create the concatenated vector */
1465: PetscCall(VecCreate(comm, &Ytmp));
1466: PetscCall(VecSetType(Ytmp, vec_type));
1467: PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1468: PetscCall(VecSetUp(Ytmp));
1469: /* copy data from X array to Y and return */
1470: for (i = 0; i < nx; i++) {
1471: PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1472: PetscCall(VecCopy(X[i], Xtmp));
1473: PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1474: }
1475: *Y = Ytmp;
1476: if (x_is) {
1477: *x_is = is_tmp;
1478: } else {
1479: for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1480: PetscCall(PetscFree(is_tmp));
1481: }
1482: }
1483: PetscFunctionReturn(PETSC_SUCCESS);
1484: }
1486: /* A helper function for VecGetSubVector to check if we can implement it with no-copy (i.e. the subvector shares
1487: memory with the original vector), and the block size of the subvector.
1489: Input Parameters:
1490: + X - the original vector
1491: - is - the index set of the subvector
1493: Output Parameters:
1494: + contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
1495: . start - start of contiguous block, as an offset from the start of the ownership range of the original vector
1496: - blocksize - the block size of the subvector
1498: */
1499: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1500: {
1501: PetscInt gstart, gend, lstart;
1502: PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1503: PetscInt n, N, ibs, vbs, bs = 1;
1505: PetscFunctionBegin;
1506: PetscCall(ISGetLocalSize(is, &n));
1507: PetscCall(ISGetSize(is, &N));
1508: PetscCall(ISGetBlockSize(is, &ibs));
1509: PetscCall(VecGetBlockSize(X, &vbs));
1510: PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1511: PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1512: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1513: if (ibs > 1) {
1514: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1515: bs = ibs;
1516: } else {
1517: if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1518: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1519: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1520: }
1522: *contig = red[0];
1523: *start = lstart;
1524: *blocksize = bs;
1525: PetscFunctionReturn(PETSC_SUCCESS);
1526: }
1528: /* A helper function for VecGetSubVector, to be used when we have to build a standalone subvector through VecScatter
1530: Input Parameters:
1531: + X - the original vector
1532: . is - the index set of the subvector
1533: - bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()
1535: Output Parameter:
1536: . Z - the subvector, which will compose the VecScatter context on output
1537: */
1538: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1539: {
1540: PetscInt n, N;
1541: VecScatter vscat;
1542: Vec Y;
1544: PetscFunctionBegin;
1545: PetscCall(ISGetLocalSize(is, &n));
1546: PetscCall(ISGetSize(is, &N));
1547: PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1548: PetscCall(VecSetSizes(Y, n, N));
1549: PetscCall(VecSetBlockSize(Y, bs));
1550: PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1551: PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1552: PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1553: PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1554: PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1555: PetscCall(VecScatterDestroy(&vscat));
1556: *Z = Y;
1557: PetscFunctionReturn(PETSC_SUCCESS);
1558: }
1560: /*@
1561: VecGetSubVector - Gets a vector representing part of another vector
1563: Collective
1565: Input Parameters:
1566: + X - vector from which to extract a subvector
1567: - is - index set representing portion of `X` to extract
1569: Output Parameter:
1570: . Y - subvector corresponding to `is`
1572: Level: advanced
1574: Notes:
1575: The subvector `Y` should be returned with `VecRestoreSubVector()`.
1576: `X` and `is` must be defined on the same communicator
1578: Changes to the subvector will be reflected in the `X` vector on the call to `VecRestoreSubVector()`.
1580: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1581: modifying the subvector. Other non-overlapping subvectors can still be obtained from `X` using this function.
1583: 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`.
1585: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1586: @*/
1587: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1588: {
1589: Vec Z;
1591: PetscFunctionBegin;
1594: PetscCheckSameComm(X, 1, is, 2);
1595: PetscAssertPointer(Y, 3);
1596: if (X->ops->getsubvector) {
1597: PetscUseTypeMethod(X, getsubvector, is, &Z);
1598: } else { /* Default implementation currently does no caching */
1599: PetscBool contig;
1600: PetscInt n, N, start, bs;
1602: PetscCall(ISGetLocalSize(is, &n));
1603: PetscCall(ISGetSize(is, &N));
1604: PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1605: if (contig) { /* We can do a no-copy implementation */
1606: const PetscScalar *x;
1607: PetscInt state = 0;
1608: PetscBool isstd, iscuda, iship;
1610: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1611: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1612: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1613: if (iscuda) {
1614: #if defined(PETSC_HAVE_CUDA)
1615: const PetscScalar *x_d;
1616: PetscMPIInt size;
1617: PetscOffloadMask flg;
1619: PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1620: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1621: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1622: if (x) x += start;
1623: if (x_d) x_d += start;
1624: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1625: if (size == 1) {
1626: PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1627: } else {
1628: PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1629: }
1630: Z->offloadmask = flg;
1631: #endif
1632: } else if (iship) {
1633: #if defined(PETSC_HAVE_HIP)
1634: const PetscScalar *x_d;
1635: PetscMPIInt size;
1636: PetscOffloadMask flg;
1638: PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1639: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1640: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1641: if (x) x += start;
1642: if (x_d) x_d += start;
1643: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1644: if (size == 1) {
1645: PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1646: } else {
1647: PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1648: }
1649: Z->offloadmask = flg;
1650: #endif
1651: } else if (isstd) {
1652: PetscMPIInt size;
1654: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1655: PetscCall(VecGetArrayRead(X, &x));
1656: if (x) x += start;
1657: if (size == 1) {
1658: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1659: } else {
1660: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1661: }
1662: PetscCall(VecRestoreArrayRead(X, &x));
1663: } else { /* default implementation: use place array */
1664: PetscCall(VecGetArrayRead(X, &x));
1665: PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1666: PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1667: PetscCall(VecSetSizes(Z, n, N));
1668: PetscCall(VecSetBlockSize(Z, bs));
1669: PetscCall(VecPlaceArray(Z, PetscSafePointerPlusOffset(x, start)));
1670: PetscCall(VecRestoreArrayRead(X, &x));
1671: }
1673: /* this is relevant only in debug mode */
1674: PetscCall(VecLockGet(X, &state));
1675: if (state) PetscCall(VecLockReadPush(Z));
1676: Z->ops->placearray = NULL;
1677: Z->ops->replacearray = NULL;
1678: } else { /* Have to create a scatter and do a copy */
1679: PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1680: }
1681: }
1682: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1683: if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1684: PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1685: *Y = Z;
1686: PetscFunctionReturn(PETSC_SUCCESS);
1687: }
1689: /*@
1690: VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`
1692: Collective
1694: Input Parameters:
1695: + X - vector from which subvector was obtained
1696: . is - index set representing the subset of `X`
1697: - Y - subvector being restored
1699: Level: advanced
1701: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1702: @*/
1703: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1704: {
1705: PETSC_UNUSED PetscObjectState dummystate = 0;
1706: PetscBool unchanged;
1708: PetscFunctionBegin;
1711: PetscCheckSameComm(X, 1, is, 2);
1712: PetscAssertPointer(Y, 3);
1715: if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1716: else {
1717: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1718: if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1719: VecScatter scatter;
1720: PetscInt state;
1722: PetscCall(VecLockGet(X, &state));
1723: PetscCheck(state == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec X is locked for read-only or read/write access");
1725: PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1726: if (scatter) {
1727: PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1728: PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1729: } else {
1730: PetscBool iscuda, iship;
1731: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1732: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1734: if (iscuda) {
1735: #if defined(PETSC_HAVE_CUDA)
1736: PetscOffloadMask ymask = (*Y)->offloadmask;
1738: /* The offloadmask of X dictates where to move memory
1739: If X GPU data is valid, then move Y data on GPU if needed
1740: Otherwise, move back to the CPU */
1741: switch (X->offloadmask) {
1742: case PETSC_OFFLOAD_BOTH:
1743: if (ymask == PETSC_OFFLOAD_CPU) {
1744: PetscCall(VecCUDAResetArray(*Y));
1745: } else if (ymask == PETSC_OFFLOAD_GPU) {
1746: X->offloadmask = PETSC_OFFLOAD_GPU;
1747: }
1748: break;
1749: case PETSC_OFFLOAD_GPU:
1750: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1751: break;
1752: case PETSC_OFFLOAD_CPU:
1753: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1754: break;
1755: case PETSC_OFFLOAD_UNALLOCATED:
1756: case PETSC_OFFLOAD_KOKKOS:
1757: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1758: }
1759: #endif
1760: } else if (iship) {
1761: #if defined(PETSC_HAVE_HIP)
1762: PetscOffloadMask ymask = (*Y)->offloadmask;
1764: /* The offloadmask of X dictates where to move memory
1765: If X GPU data is valid, then move Y data on GPU if needed
1766: Otherwise, move back to the CPU */
1767: switch (X->offloadmask) {
1768: case PETSC_OFFLOAD_BOTH:
1769: if (ymask == PETSC_OFFLOAD_CPU) {
1770: PetscCall(VecHIPResetArray(*Y));
1771: } else if (ymask == PETSC_OFFLOAD_GPU) {
1772: X->offloadmask = PETSC_OFFLOAD_GPU;
1773: }
1774: break;
1775: case PETSC_OFFLOAD_GPU:
1776: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1777: break;
1778: case PETSC_OFFLOAD_CPU:
1779: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1780: break;
1781: case PETSC_OFFLOAD_UNALLOCATED:
1782: case PETSC_OFFLOAD_KOKKOS:
1783: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1784: }
1785: #endif
1786: } else {
1787: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1788: PetscCall(VecResetArray(*Y));
1789: }
1790: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1791: }
1792: }
1793: }
1794: PetscCall(VecDestroy(Y));
1795: PetscFunctionReturn(PETSC_SUCCESS);
1796: }
1798: /*@
1799: VecCreateLocalVector - Creates a vector object suitable for use with `VecGetLocalVector()` and friends. You must call `VecDestroy()` when the
1800: vector is no longer needed.
1802: Not Collective.
1804: Input Parameter:
1805: . v - The vector for which the local vector is desired.
1807: Output Parameter:
1808: . w - Upon exit this contains the local vector.
1810: Level: beginner
1812: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1813: @*/
1814: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1815: {
1816: VecType roottype;
1817: PetscInt n;
1819: PetscFunctionBegin;
1821: PetscAssertPointer(w, 2);
1822: if (v->ops->createlocalvector) {
1823: PetscUseTypeMethod(v, createlocalvector, w);
1824: PetscFunctionReturn(PETSC_SUCCESS);
1825: }
1826: PetscCall(VecGetRootType_Private(v, &roottype));
1827: PetscCall(VecCreate(PETSC_COMM_SELF, w));
1828: PetscCall(VecGetLocalSize(v, &n));
1829: PetscCall(VecSetSizes(*w, n, n));
1830: PetscCall(VecGetBlockSize(v, &n));
1831: PetscCall(VecSetBlockSize(*w, n));
1832: PetscCall(VecSetType(*w, roottype));
1833: PetscFunctionReturn(PETSC_SUCCESS);
1834: }
1836: /*@
1837: VecGetLocalVectorRead - Maps the local portion of a vector into a
1838: vector.
1840: Not Collective.
1842: Input Parameter:
1843: . v - The vector for which the local vector is desired.
1845: Output Parameter:
1846: . w - Upon exit this contains the local vector.
1848: Level: beginner
1850: Notes:
1851: You must call `VecRestoreLocalVectorRead()` when the local
1852: vector is no longer needed.
1854: This function is similar to `VecGetArrayRead()` which maps the local
1855: portion into a raw pointer. `VecGetLocalVectorRead()` is usually
1856: almost as efficient as `VecGetArrayRead()` but in certain circumstances
1857: `VecGetLocalVectorRead()` can be much more efficient than
1858: `VecGetArrayRead()`. This is because the construction of a contiguous
1859: array representing the vector data required by `VecGetArrayRead()` can
1860: be an expensive operation for certain vector types. For example, for
1861: GPU vectors `VecGetArrayRead()` requires that the data between device
1862: and host is synchronized.
1864: Unlike `VecGetLocalVector()`, this routine is not collective and
1865: preserves cached information.
1867: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1868: @*/
1869: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1870: {
1871: PetscFunctionBegin;
1874: VecCheckSameLocalSize(v, 1, w, 2);
1875: if (v->ops->getlocalvectorread) {
1876: PetscUseTypeMethod(v, getlocalvectorread, w);
1877: } else {
1878: PetscScalar *a;
1880: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1881: PetscCall(VecPlaceArray(w, a));
1882: }
1883: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1884: PetscCall(VecLockReadPush(v));
1885: PetscCall(VecLockReadPush(w));
1886: PetscFunctionReturn(PETSC_SUCCESS);
1887: }
1889: /*@
1890: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1891: previously mapped into a vector using `VecGetLocalVectorRead()`.
1893: Not Collective.
1895: Input Parameters:
1896: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
1897: - w - The vector into which the local portion of `v` was mapped.
1899: Level: beginner
1901: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1902: @*/
1903: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1904: {
1905: PetscFunctionBegin;
1908: if (v->ops->restorelocalvectorread) {
1909: PetscUseTypeMethod(v, restorelocalvectorread, w);
1910: } else {
1911: const PetscScalar *a;
1913: PetscCall(VecGetArrayRead(w, &a));
1914: PetscCall(VecRestoreArrayRead(v, &a));
1915: PetscCall(VecResetArray(w));
1916: }
1917: PetscCall(VecLockReadPop(v));
1918: PetscCall(VecLockReadPop(w));
1919: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1920: PetscFunctionReturn(PETSC_SUCCESS);
1921: }
1923: /*@
1924: VecGetLocalVector - Maps the local portion of a vector into a
1925: vector.
1927: Collective
1929: Input Parameter:
1930: . v - The vector for which the local vector is desired.
1932: Output Parameter:
1933: . w - Upon exit this contains the local vector.
1935: Level: beginner
1937: Notes:
1938: You must call `VecRestoreLocalVector()` when the local
1939: vector is no longer needed.
1941: This function is similar to `VecGetArray()` which maps the local
1942: portion into a raw pointer. `VecGetLocalVector()` is usually about as
1943: efficient as `VecGetArray()` but in certain circumstances
1944: `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1945: This is because the construction of a contiguous array representing
1946: the vector data required by `VecGetArray()` can be an expensive
1947: operation for certain vector types. For example, for GPU vectors
1948: `VecGetArray()` requires that the data between device and host is
1949: synchronized.
1951: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1952: @*/
1953: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1954: {
1955: PetscFunctionBegin;
1958: VecCheckSameLocalSize(v, 1, w, 2);
1959: if (v->ops->getlocalvector) {
1960: PetscUseTypeMethod(v, getlocalvector, w);
1961: } else {
1962: PetscScalar *a;
1964: PetscCall(VecGetArray(v, &a));
1965: PetscCall(VecPlaceArray(w, a));
1966: }
1967: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1968: PetscFunctionReturn(PETSC_SUCCESS);
1969: }
1971: /*@
1972: VecRestoreLocalVector - Unmaps the local portion of a vector
1973: previously mapped into a vector using `VecGetLocalVector()`.
1975: Logically Collective.
1977: Input Parameters:
1978: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
1979: - w - The vector into which the local portion of `v` was mapped.
1981: Level: beginner
1983: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1984: @*/
1985: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1986: {
1987: PetscFunctionBegin;
1990: if (v->ops->restorelocalvector) {
1991: PetscUseTypeMethod(v, restorelocalvector, w);
1992: } else {
1993: PetscScalar *a;
1994: PetscCall(VecGetArray(w, &a));
1995: PetscCall(VecRestoreArray(v, &a));
1996: PetscCall(VecResetArray(w));
1997: }
1998: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1999: PetscCall(PetscObjectStateIncrease((PetscObject)v));
2000: PetscFunctionReturn(PETSC_SUCCESS);
2001: }
2003: /*@C
2004: VecGetArray - Returns a pointer to a contiguous array that contains this
2005: MPI processes's portion of the vector data
2007: Logically Collective
2009: Input Parameter:
2010: . x - the vector
2012: Output Parameter:
2013: . a - location to put pointer to the array
2015: Level: beginner
2017: Notes:
2018: For the standard PETSc vectors, `VecGetArray()` returns a pointer to the local data array and
2019: does not use any copies. If the underlying vector data is not stored in a contiguous array
2020: this routine will copy the data to a contiguous array and return a pointer to that. You MUST
2021: call `VecRestoreArray()` when you no longer need access to the array.
2023: Fortran Note:
2024: .vb
2025: PetscScalar, pointer :: a(:)
2026: .ve
2028: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecPlaceArray()`, `VecGetArray2d()`,
2029: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2030: @*/
2031: PetscErrorCode VecGetArray(Vec x, PetscScalar *a[])
2032: {
2033: PetscFunctionBegin;
2035: PetscCall(VecSetErrorIfLocked(x, 1));
2036: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
2037: PetscUseTypeMethod(x, getarray, a);
2038: } else if (x->petscnative) { /* VECSTANDARD */
2039: *a = *((PetscScalar **)x->data);
2040: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
2041: PetscFunctionReturn(PETSC_SUCCESS);
2042: }
2044: /*@C
2045: VecRestoreArray - Restores a vector after `VecGetArray()` has been called and the array is no longer needed
2047: Logically Collective
2049: Input Parameters:
2050: + x - the vector
2051: - a - location of pointer to array obtained from `VecGetArray()`
2053: Level: beginner
2055: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2056: `VecGetArrayPair()`, `VecRestoreArrayPair()`
2057: @*/
2058: PetscErrorCode VecRestoreArray(Vec x, PetscScalar *a[])
2059: {
2060: PetscFunctionBegin;
2062: if (a) PetscAssertPointer(a, 2);
2063: if (x->ops->restorearray) {
2064: PetscUseTypeMethod(x, restorearray, a);
2065: } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
2066: if (a) *a = NULL;
2067: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2068: PetscFunctionReturn(PETSC_SUCCESS);
2069: }
2070: /*@C
2071: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
2073: Not Collective
2075: Input Parameter:
2076: . x - the vector
2078: Output Parameter:
2079: . a - the array
2081: Level: beginner
2083: Notes:
2084: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2086: Unlike `VecGetArray()`, preserves cached information like vector norms.
2088: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
2089: implementations may require a copy, but such implementations should cache the contiguous representation so that
2090: only one copy is performed when this routine is called multiple times in sequence.
2092: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2093: @*/
2094: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar *a[])
2095: {
2096: PetscFunctionBegin;
2098: PetscAssertPointer(a, 2);
2099: if (x->ops->getarrayread) {
2100: PetscUseTypeMethod(x, getarrayread, a);
2101: } else if (x->ops->getarray) {
2102: PetscObjectState state;
2104: /* VECNEST, VECCUDA, VECKOKKOS etc */
2105: // x->ops->getarray may bump the object state, but since we know this is a read-only get
2106: // we can just undo that
2107: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2108: PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
2109: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2110: } else if (x->petscnative) {
2111: /* VECSTANDARD */
2112: *a = *((PetscScalar **)x->data);
2113: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
2114: PetscFunctionReturn(PETSC_SUCCESS);
2115: }
2117: /*@C
2118: VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`
2120: Not Collective
2122: Input Parameters:
2123: + x - the vector
2124: - a - the array
2126: Level: beginner
2128: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2129: @*/
2130: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar *a[])
2131: {
2132: PetscFunctionBegin;
2134: if (a) PetscAssertPointer(a, 2);
2135: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
2136: /* nothing */
2137: } else if (x->ops->restorearrayread) { /* VECNEST */
2138: PetscUseTypeMethod(x, restorearrayread, a);
2139: } else { /* No one? */
2140: PetscObjectState state;
2142: // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2143: // we can just undo that
2144: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2145: PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2146: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2147: }
2148: if (a) *a = NULL;
2149: PetscFunctionReturn(PETSC_SUCCESS);
2150: }
2152: /*@C
2153: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
2154: MPI processes's portion of the vector data.
2156: Logically Collective
2158: Input Parameter:
2159: . x - the vector
2161: Output Parameter:
2162: . a - location to put pointer to the array
2164: Level: intermediate
2166: Note:
2167: The values in this array are NOT valid, the caller of this routine is responsible for putting
2168: values into the array; any values it does not set will be invalid.
2170: The array must be returned using a matching call to `VecRestoreArrayWrite()`.
2172: For vectors associated with GPUs, the host and device vectors are not synchronized before
2173: giving access. If you need correct values in the array use `VecGetArray()`
2175: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecPlaceArray()`, `VecGetArray2d()`,
2176: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
2177: @*/
2178: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar *a[])
2179: {
2180: PetscFunctionBegin;
2182: PetscAssertPointer(a, 2);
2183: PetscCall(VecSetErrorIfLocked(x, 1));
2184: if (x->ops->getarraywrite) {
2185: PetscUseTypeMethod(x, getarraywrite, a);
2186: } else {
2187: PetscCall(VecGetArray(x, a));
2188: }
2189: PetscFunctionReturn(PETSC_SUCCESS);
2190: }
2192: /*@C
2193: VecRestoreArrayWrite - Restores a vector after `VecGetArrayWrite()` has been called.
2195: Logically Collective
2197: Input Parameters:
2198: + x - the vector
2199: - a - location of pointer to array obtained from `VecGetArray()`
2201: Level: beginner
2203: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2204: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2205: @*/
2206: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar *a[])
2207: {
2208: PetscFunctionBegin;
2210: if (a) PetscAssertPointer(a, 2);
2211: if (x->ops->restorearraywrite) {
2212: PetscUseTypeMethod(x, restorearraywrite, a);
2213: } else if (x->ops->restorearray) {
2214: PetscUseTypeMethod(x, restorearray, a);
2215: }
2216: if (a) *a = NULL;
2217: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2218: PetscFunctionReturn(PETSC_SUCCESS);
2219: }
2221: /*@C
2222: VecGetArrays - Returns a pointer to the arrays in a set of vectors
2223: that were created by a call to `VecDuplicateVecs()`.
2225: Logically Collective; No Fortran Support
2227: Input Parameters:
2228: + x - the vectors
2229: - n - the number of vectors
2231: Output Parameter:
2232: . a - location to put pointer to the array
2234: Level: intermediate
2236: Note:
2237: You MUST call `VecRestoreArrays()` when you no longer need access to the arrays.
2239: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2240: @*/
2241: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2242: {
2243: PetscInt i;
2244: PetscScalar **q;
2246: PetscFunctionBegin;
2247: PetscAssertPointer(x, 1);
2249: PetscAssertPointer(a, 3);
2250: PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2251: PetscCall(PetscMalloc1(n, &q));
2252: for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2253: *a = q;
2254: PetscFunctionReturn(PETSC_SUCCESS);
2255: }
2257: /*@C
2258: VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2259: has been called.
2261: Logically Collective; No Fortran Support
2263: Input Parameters:
2264: + x - the vector
2265: . n - the number of vectors
2266: - a - location of pointer to arrays obtained from `VecGetArrays()`
2268: Notes:
2269: For regular PETSc vectors this routine does not involve any copies. For
2270: any special vectors that do not store local vector data in a contiguous
2271: array, this routine will copy the data back into the underlying
2272: vector data structure from the arrays obtained with `VecGetArrays()`.
2274: Level: intermediate
2276: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2277: @*/
2278: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2279: {
2280: PetscInt i;
2281: PetscScalar **q = *a;
2283: PetscFunctionBegin;
2284: PetscAssertPointer(x, 1);
2286: PetscAssertPointer(a, 3);
2288: for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2289: PetscCall(PetscFree(q));
2290: PetscFunctionReturn(PETSC_SUCCESS);
2291: }
2293: /*@C
2294: VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g.,
2295: `VECCUDA`), the returned pointer will be a device pointer to the device memory that contains
2296: this MPI processes's portion of the vector data.
2298: Logically Collective; No Fortran Support
2300: Input Parameter:
2301: . x - the vector
2303: Output Parameters:
2304: + a - location to put pointer to the array
2305: - mtype - memory type of the array
2307: Level: beginner
2309: Note:
2310: Device data is guaranteed to have the latest value. Otherwise, when this is a host vector
2311: (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host
2312: pointer.
2314: For `VECKOKKOS`, if Kokkos is configured without device (e.g., use serial or openmp), per
2315: this function, the vector works like `VECSEQ`/`VECMPI`; otherwise, it works like `VECCUDA` or
2316: `VECHIP` etc.
2318: Use `VecRestoreArrayAndMemType()` when the array access is no longer needed.
2320: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`,
2321: `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2322: @*/
2323: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar *a[], PetscMemType *mtype)
2324: {
2325: PetscFunctionBegin;
2328: PetscAssertPointer(a, 2);
2329: if (mtype) PetscAssertPointer(mtype, 3);
2330: PetscCall(VecSetErrorIfLocked(x, 1));
2331: if (x->ops->getarrayandmemtype) {
2332: /* VECCUDA, VECKOKKOS etc */
2333: PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2334: } else {
2335: /* VECSTANDARD, VECNEST, VECVIENNACL */
2336: PetscCall(VecGetArray(x, a));
2337: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2338: }
2339: PetscFunctionReturn(PETSC_SUCCESS);
2340: }
2342: /*@C
2343: VecRestoreArrayAndMemType - Restores a vector after `VecGetArrayAndMemType()` has been called.
2345: Logically Collective; No Fortran Support
2347: Input Parameters:
2348: + x - the vector
2349: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`
2351: Level: beginner
2353: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`,
2354: `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2355: @*/
2356: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar *a[])
2357: {
2358: PetscFunctionBegin;
2361: if (a) PetscAssertPointer(a, 2);
2362: if (x->ops->restorearrayandmemtype) {
2363: /* VECCUDA, VECKOKKOS etc */
2364: PetscUseTypeMethod(x, restorearrayandmemtype, a);
2365: } else {
2366: /* VECNEST, VECVIENNACL */
2367: PetscCall(VecRestoreArray(x, a));
2368: } /* VECSTANDARD does nothing */
2369: if (a) *a = NULL;
2370: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2371: PetscFunctionReturn(PETSC_SUCCESS);
2372: }
2374: /*@C
2375: VecGetArrayReadAndMemType - Like `VecGetArrayRead()`, but if the input vector is a device vector, it will return a read-only device pointer.
2376: The returned pointer is guaranteed to point to up-to-date data. For host vectors, it functions as `VecGetArrayRead()`.
2378: Not Collective; No Fortran Support
2380: Input Parameter:
2381: . x - the vector
2383: Output Parameters:
2384: + a - the array
2385: - mtype - memory type of the array
2387: Level: beginner
2389: Notes:
2390: The array must be returned using a matching call to `VecRestoreArrayReadAndMemType()`.
2392: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2393: @*/
2394: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar *a[], PetscMemType *mtype)
2395: {
2396: PetscFunctionBegin;
2399: PetscAssertPointer(a, 2);
2400: if (mtype) PetscAssertPointer(mtype, 3);
2401: if (x->ops->getarrayreadandmemtype) {
2402: /* VECCUDA/VECHIP though they are also petscnative */
2403: PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2404: } else if (x->ops->getarrayandmemtype) {
2405: /* VECKOKKOS */
2406: PetscObjectState state;
2408: // see VecGetArrayRead() for why
2409: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2410: PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2411: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2412: } else {
2413: PetscCall(VecGetArrayRead(x, a));
2414: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2415: }
2416: PetscFunctionReturn(PETSC_SUCCESS);
2417: }
2419: /*@C
2420: VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`
2422: Not Collective; No Fortran Support
2424: Input Parameters:
2425: + x - the vector
2426: - a - the array
2428: Level: beginner
2430: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2431: @*/
2432: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar *a[])
2433: {
2434: PetscFunctionBegin;
2437: if (a) PetscAssertPointer(a, 2);
2438: if (x->ops->restorearrayreadandmemtype) {
2439: /* VECCUDA/VECHIP */
2440: PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2441: } else if (!x->petscnative) {
2442: /* VECNEST */
2443: PetscCall(VecRestoreArrayRead(x, a));
2444: }
2445: if (a) *a = NULL;
2446: PetscFunctionReturn(PETSC_SUCCESS);
2447: }
2449: /*@C
2450: VecGetArrayWriteAndMemType - Like `VecGetArrayWrite()`, but if this is a device vector it will always return
2451: a device pointer to the device memory that contains this processor's portion of the vector data.
2453: Logically Collective; No Fortran Support
2455: Input Parameter:
2456: . x - the vector
2458: Output Parameters:
2459: + a - the array
2460: - mtype - memory type of the array
2462: Level: beginner
2464: Note:
2465: The array must be returned using a matching call to `VecRestoreArrayWriteAndMemType()`, where it will label the device memory as most recent.
2467: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2468: @*/
2469: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar *a[], PetscMemType *mtype)
2470: {
2471: PetscFunctionBegin;
2474: PetscCall(VecSetErrorIfLocked(x, 1));
2475: PetscAssertPointer(a, 2);
2476: if (mtype) PetscAssertPointer(mtype, 3);
2477: if (x->ops->getarraywriteandmemtype) {
2478: /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2479: PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2480: } else if (x->ops->getarrayandmemtype) {
2481: PetscCall(VecGetArrayAndMemType(x, a, mtype));
2482: } else {
2483: /* VECNEST, VECVIENNACL */
2484: PetscCall(VecGetArrayWrite(x, a));
2485: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2486: }
2487: PetscFunctionReturn(PETSC_SUCCESS);
2488: }
2490: /*@C
2491: VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`
2493: Logically Collective; No Fortran Support
2495: Input Parameters:
2496: + x - the vector
2497: - a - the array
2499: Level: beginner
2501: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2502: @*/
2503: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar *a[])
2504: {
2505: PetscFunctionBegin;
2508: PetscCall(VecSetErrorIfLocked(x, 1));
2509: if (a) PetscAssertPointer(a, 2);
2510: if (x->ops->restorearraywriteandmemtype) {
2511: /* VECCUDA/VECHIP */
2512: PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2513: PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2514: } else if (x->ops->restorearrayandmemtype) {
2515: PetscCall(VecRestoreArrayAndMemType(x, a));
2516: } else {
2517: PetscCall(VecRestoreArray(x, a));
2518: }
2519: if (a) *a = NULL;
2520: PetscFunctionReturn(PETSC_SUCCESS);
2521: }
2523: /*@
2524: VecPlaceArray - Allows one to replace the array in a vector with an
2525: array provided by the user. This is useful to avoid copying an array
2526: into a vector.
2528: Logically Collective; No Fortran Support
2530: Input Parameters:
2531: + vec - the vector
2532: - array - the array
2534: Level: developer
2536: Notes:
2537: Adding `const` to `array` was an oversight, as subsequent operations on `vec` would
2538: likely modify the data in `array`. However, we have kept it to avoid breaking APIs.
2540: Use `VecReplaceArray()` instead to permanently replace the array
2542: You can return to the original array with a call to `VecResetArray()`. `vec` does not take
2543: ownership of `array` in any way.
2545: The user must free `array` themselves but be careful not to
2546: do so before the vector has either been destroyed, had its original array restored with
2547: `VecResetArray()` or permanently replaced with `VecReplaceArray()`.
2549: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2550: @*/
2551: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2552: {
2553: PetscFunctionBegin;
2556: if (array) PetscAssertPointer(array, 2);
2557: PetscUseTypeMethod(vec, placearray, array);
2558: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2559: PetscFunctionReturn(PETSC_SUCCESS);
2560: }
2562: /*@C
2563: VecReplaceArray - Allows one to replace the array in a vector with an
2564: array provided by the user. This is useful to avoid copying an array
2565: into a vector.
2567: Logically Collective; No Fortran Support
2569: Input Parameters:
2570: + vec - the vector
2571: - array - the array
2573: Level: developer
2575: Notes:
2576: Adding `const` to `array` was an oversight, as subsequent operations on `vec` would
2577: likely modify the data in `array`. However, we have kept it to avoid breaking APIs.
2579: This permanently replaces the array and frees the memory associated
2580: with the old array. Use `VecPlaceArray()` to temporarily replace the array.
2582: The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
2583: freed by the user. It will be freed when the vector is destroyed.
2585: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2586: @*/
2587: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2588: {
2589: PetscFunctionBegin;
2592: PetscUseTypeMethod(vec, replacearray, array);
2593: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2594: PetscFunctionReturn(PETSC_SUCCESS);
2595: }
2597: /*@C
2598: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2599: processor's portion of the vector data. You MUST call `VecRestoreArray2d()`
2600: when you no longer need access to the array.
2602: Logically Collective
2604: Input Parameters:
2605: + x - the vector
2606: . m - first dimension of two dimensional array
2607: . n - second dimension of two dimensional array
2608: . mstart - first index you will use in first coordinate direction (often 0)
2609: - nstart - first index in the second coordinate direction (often 0)
2611: Output Parameter:
2612: . a - location to put pointer to the array
2614: Level: developer
2616: Notes:
2617: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2618: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2619: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2620: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2622: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2624: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2625: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2626: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2627: @*/
2628: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2629: {
2630: PetscInt i, N;
2631: PetscScalar *aa;
2633: PetscFunctionBegin;
2635: PetscAssertPointer(a, 6);
2637: PetscCall(VecGetLocalSize(x, &N));
2638: 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);
2639: PetscCall(VecGetArray(x, &aa));
2641: PetscCall(PetscMalloc1(m, a));
2642: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2643: *a -= mstart;
2644: PetscFunctionReturn(PETSC_SUCCESS);
2645: }
2647: /*@C
2648: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2649: processor's portion of the vector data. You MUST call `VecRestoreArray2dWrite()`
2650: when you no longer need access to the array.
2652: Logically Collective
2654: Input Parameters:
2655: + x - the vector
2656: . m - first dimension of two dimensional array
2657: . n - second dimension of two dimensional array
2658: . mstart - first index you will use in first coordinate direction (often 0)
2659: - nstart - first index in the second coordinate direction (often 0)
2661: Output Parameter:
2662: . a - location to put pointer to the array
2664: Level: developer
2666: Notes:
2667: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2668: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2669: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2670: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2672: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2674: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2675: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2676: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2677: @*/
2678: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2679: {
2680: PetscInt i, N;
2681: PetscScalar *aa;
2683: PetscFunctionBegin;
2685: PetscAssertPointer(a, 6);
2687: PetscCall(VecGetLocalSize(x, &N));
2688: 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);
2689: PetscCall(VecGetArrayWrite(x, &aa));
2691: PetscCall(PetscMalloc1(m, a));
2692: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2693: *a -= mstart;
2694: PetscFunctionReturn(PETSC_SUCCESS);
2695: }
2697: /*@C
2698: VecRestoreArray2d - Restores a vector after `VecGetArray2d()` has been called.
2700: Logically Collective
2702: Input Parameters:
2703: + x - the vector
2704: . m - first dimension of two dimensional array
2705: . n - second dimension of the two dimensional array
2706: . mstart - first index you will use in first coordinate direction (often 0)
2707: . nstart - first index in the second coordinate direction (often 0)
2708: - a - location of pointer to array obtained from `VecGetArray2d()`
2710: Level: developer
2712: Notes:
2713: For regular PETSc vectors this routine does not involve any copies. For
2714: any special vectors that do not store local vector data in a contiguous
2715: array, this routine will copy the data back into the underlying
2716: vector data structure from the array obtained with `VecGetArray()`.
2718: This routine actually zeros out the `a` pointer.
2720: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2721: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2722: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2723: @*/
2724: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2725: {
2726: void *dummy;
2728: PetscFunctionBegin;
2730: PetscAssertPointer(a, 6);
2732: dummy = (void *)(*a + mstart);
2733: PetscCall(PetscFree(dummy));
2734: PetscCall(VecRestoreArray(x, NULL));
2735: *a = NULL;
2736: PetscFunctionReturn(PETSC_SUCCESS);
2737: }
2739: /*@C
2740: VecRestoreArray2dWrite - Restores a vector after `VecGetArray2dWrite()` has been called.
2742: Logically Collective
2744: Input Parameters:
2745: + x - the vector
2746: . m - first dimension of two dimensional array
2747: . n - second dimension of the two dimensional array
2748: . mstart - first index you will use in first coordinate direction (often 0)
2749: . nstart - first index in the second coordinate direction (often 0)
2750: - a - location of pointer to array obtained from `VecGetArray2d()`
2752: Level: developer
2754: Notes:
2755: For regular PETSc vectors this routine does not involve any copies. For
2756: any special vectors that do not store local vector data in a contiguous
2757: array, this routine will copy the data back into the underlying
2758: vector data structure from the array obtained with `VecGetArray()`.
2760: This routine actually zeros out the `a` pointer.
2762: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2763: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2764: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2765: @*/
2766: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2767: {
2768: void *dummy;
2770: PetscFunctionBegin;
2772: PetscAssertPointer(a, 6);
2774: dummy = (void *)(*a + mstart);
2775: PetscCall(PetscFree(dummy));
2776: PetscCall(VecRestoreArrayWrite(x, NULL));
2777: PetscFunctionReturn(PETSC_SUCCESS);
2778: }
2780: /*@C
2781: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2782: processor's portion of the vector data. You MUST call `VecRestoreArray1d()`
2783: when you no longer need access to the array.
2785: Logically Collective
2787: Input Parameters:
2788: + x - the vector
2789: . m - first dimension of two dimensional array
2790: - mstart - first index you will use in first coordinate direction (often 0)
2792: Output Parameter:
2793: . a - location to put pointer to the array
2795: Level: developer
2797: Notes:
2798: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2799: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2800: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2802: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2804: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2805: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2806: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2807: @*/
2808: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2809: {
2810: PetscInt N;
2812: PetscFunctionBegin;
2814: PetscAssertPointer(a, 4);
2816: PetscCall(VecGetLocalSize(x, &N));
2817: 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);
2818: PetscCall(VecGetArray(x, a));
2819: *a -= mstart;
2820: PetscFunctionReturn(PETSC_SUCCESS);
2821: }
2823: /*@C
2824: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
2825: processor's portion of the vector data. You MUST call `VecRestoreArray1dWrite()`
2826: when you no longer need access to the array.
2828: Logically Collective
2830: Input Parameters:
2831: + x - the vector
2832: . m - first dimension of two dimensional array
2833: - mstart - first index you will use in first coordinate direction (often 0)
2835: Output Parameter:
2836: . a - location to put pointer to the array
2838: Level: developer
2840: Notes:
2841: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2842: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2843: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2845: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2847: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2848: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2849: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2850: @*/
2851: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2852: {
2853: PetscInt N;
2855: PetscFunctionBegin;
2857: PetscAssertPointer(a, 4);
2859: PetscCall(VecGetLocalSize(x, &N));
2860: 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);
2861: PetscCall(VecGetArrayWrite(x, a));
2862: *a -= mstart;
2863: PetscFunctionReturn(PETSC_SUCCESS);
2864: }
2866: /*@C
2867: VecRestoreArray1d - Restores a vector after `VecGetArray1d()` has been called.
2869: Logically Collective
2871: Input Parameters:
2872: + x - the vector
2873: . m - first dimension of two dimensional array
2874: . mstart - first index you will use in first coordinate direction (often 0)
2875: - a - location of pointer to array obtained from `VecGetArray1d()`
2877: Level: developer
2879: Notes:
2880: For regular PETSc vectors this routine does not involve any copies. For
2881: any special vectors that do not store local vector data in a contiguous
2882: array, this routine will copy the data back into the underlying
2883: vector data structure from the array obtained with `VecGetArray1d()`.
2885: This routine actually zeros out the `a` pointer.
2887: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2888: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2889: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2890: @*/
2891: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2892: {
2893: PetscFunctionBegin;
2896: PetscCall(VecRestoreArray(x, NULL));
2897: *a = NULL;
2898: PetscFunctionReturn(PETSC_SUCCESS);
2899: }
2901: /*@C
2902: VecRestoreArray1dWrite - Restores a vector after `VecGetArray1dWrite()` has been called.
2904: Logically Collective
2906: Input Parameters:
2907: + x - the vector
2908: . m - first dimension of two dimensional array
2909: . mstart - first index you will use in first coordinate direction (often 0)
2910: - a - location of pointer to array obtained from `VecGetArray1d()`
2912: Level: developer
2914: Notes:
2915: For regular PETSc vectors this routine does not involve any copies. For
2916: any special vectors that do not store local vector data in a contiguous
2917: array, this routine will copy the data back into the underlying
2918: vector data structure from the array obtained with `VecGetArray1d()`.
2920: This routine actually zeros out the `a` pointer.
2922: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2923: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2924: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2925: @*/
2926: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2927: {
2928: PetscFunctionBegin;
2931: PetscCall(VecRestoreArrayWrite(x, NULL));
2932: *a = NULL;
2933: PetscFunctionReturn(PETSC_SUCCESS);
2934: }
2936: /*@C
2937: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
2938: processor's portion of the vector data. You MUST call `VecRestoreArray3d()`
2939: when you no longer need access to the array.
2941: Logically Collective
2943: Input Parameters:
2944: + x - the vector
2945: . m - first dimension of three dimensional array
2946: . n - second dimension of three dimensional array
2947: . p - third dimension of three dimensional array
2948: . mstart - first index you will use in first coordinate direction (often 0)
2949: . nstart - first index in the second coordinate direction (often 0)
2950: - pstart - first index in the third coordinate direction (often 0)
2952: Output Parameter:
2953: . a - location to put pointer to the array
2955: Level: developer
2957: Notes:
2958: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
2959: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2960: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2961: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
2963: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2965: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2966: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
2967: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2968: @*/
2969: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
2970: {
2971: PetscInt i, N, j;
2972: PetscScalar *aa, **b;
2974: PetscFunctionBegin;
2976: PetscAssertPointer(a, 8);
2978: PetscCall(VecGetLocalSize(x, &N));
2979: 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);
2980: PetscCall(VecGetArray(x, &aa));
2982: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
2983: b = (PetscScalar **)((*a) + m);
2984: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
2985: for (i = 0; i < m; i++)
2986: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset(aa, i * n * p + j * p - pstart);
2987: *a -= mstart;
2988: PetscFunctionReturn(PETSC_SUCCESS);
2989: }
2991: /*@C
2992: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
2993: processor's portion of the vector data. You MUST call `VecRestoreArray3dWrite()`
2994: when you no longer need access to the array.
2996: Logically Collective
2998: Input Parameters:
2999: + x - the vector
3000: . m - first dimension of three dimensional array
3001: . n - second dimension of three dimensional array
3002: . p - third dimension of three dimensional array
3003: . mstart - first index you will use in first coordinate direction (often 0)
3004: . nstart - first index in the second coordinate direction (often 0)
3005: - pstart - first index in the third coordinate direction (often 0)
3007: Output Parameter:
3008: . a - location to put pointer to the array
3010: Level: developer
3012: Notes:
3013: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3014: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3015: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3016: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3018: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3020: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3021: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3022: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3023: @*/
3024: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3025: {
3026: PetscInt i, N, j;
3027: PetscScalar *aa, **b;
3029: PetscFunctionBegin;
3031: PetscAssertPointer(a, 8);
3033: PetscCall(VecGetLocalSize(x, &N));
3034: 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);
3035: PetscCall(VecGetArrayWrite(x, &aa));
3037: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3038: b = (PetscScalar **)((*a) + m);
3039: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3040: for (i = 0; i < m; i++)
3041: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3043: *a -= mstart;
3044: PetscFunctionReturn(PETSC_SUCCESS);
3045: }
3047: /*@C
3048: VecRestoreArray3d - Restores a vector after `VecGetArray3d()` has been called.
3050: Logically Collective
3052: Input Parameters:
3053: + x - the vector
3054: . m - first dimension of three dimensional array
3055: . n - second dimension of the three dimensional array
3056: . p - third dimension of the three dimensional array
3057: . mstart - first index you will use in first coordinate direction (often 0)
3058: . nstart - first index in the second coordinate direction (often 0)
3059: . pstart - first index in the third coordinate direction (often 0)
3060: - a - location of pointer to array obtained from VecGetArray3d()
3062: Level: developer
3064: Notes:
3065: For regular PETSc vectors this routine does not involve any copies. For
3066: any special vectors that do not store local vector data in a contiguous
3067: array, this routine will copy the data back into the underlying
3068: vector data structure from the array obtained with `VecGetArray()`.
3070: This routine actually zeros out the `a` pointer.
3072: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3073: `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3074: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3075: @*/
3076: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3077: {
3078: void *dummy;
3080: PetscFunctionBegin;
3082: PetscAssertPointer(a, 8);
3084: dummy = (void *)(*a + mstart);
3085: PetscCall(PetscFree(dummy));
3086: PetscCall(VecRestoreArray(x, NULL));
3087: *a = NULL;
3088: PetscFunctionReturn(PETSC_SUCCESS);
3089: }
3091: /*@C
3092: VecRestoreArray3dWrite - Restores a vector after `VecGetArray3dWrite()` has been called.
3094: Logically Collective
3096: Input Parameters:
3097: + x - the vector
3098: . m - first dimension of three dimensional array
3099: . n - second dimension of the three dimensional array
3100: . p - third dimension of the three dimensional array
3101: . mstart - first index you will use in first coordinate direction (often 0)
3102: . nstart - first index in the second coordinate direction (often 0)
3103: . pstart - first index in the third coordinate direction (often 0)
3104: - a - location of pointer to array obtained from VecGetArray3d()
3106: Level: developer
3108: Notes:
3109: For regular PETSc vectors this routine does not involve any copies. For
3110: any special vectors that do not store local vector data in a contiguous
3111: array, this routine will copy the data back into the underlying
3112: vector data structure from the array obtained with `VecGetArray()`.
3114: This routine actually zeros out the `a` pointer.
3116: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3117: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3118: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3119: @*/
3120: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3121: {
3122: void *dummy;
3124: PetscFunctionBegin;
3126: PetscAssertPointer(a, 8);
3128: dummy = (void *)(*a + mstart);
3129: PetscCall(PetscFree(dummy));
3130: PetscCall(VecRestoreArrayWrite(x, NULL));
3131: *a = NULL;
3132: PetscFunctionReturn(PETSC_SUCCESS);
3133: }
3135: /*@C
3136: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3137: processor's portion of the vector data. You MUST call `VecRestoreArray4d()`
3138: when you no longer need access to the array.
3140: Logically Collective
3142: Input Parameters:
3143: + x - the vector
3144: . m - first dimension of four dimensional array
3145: . n - second dimension of four dimensional array
3146: . p - third dimension of four dimensional array
3147: . q - fourth dimension of four dimensional array
3148: . mstart - first index you will use in first coordinate direction (often 0)
3149: . nstart - first index in the second coordinate direction (often 0)
3150: . pstart - first index in the third coordinate direction (often 0)
3151: - qstart - first index in the fourth coordinate direction (often 0)
3153: Output Parameter:
3154: . a - location to put pointer to the array
3156: Level: developer
3158: Notes:
3159: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3160: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3161: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3162: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3164: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3166: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3167: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3168: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3169: @*/
3170: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3171: {
3172: PetscInt i, N, j, k;
3173: PetscScalar *aa, ***b, **c;
3175: PetscFunctionBegin;
3177: PetscAssertPointer(a, 10);
3179: PetscCall(VecGetLocalSize(x, &N));
3180: 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);
3181: PetscCall(VecGetArray(x, &aa));
3183: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3184: b = (PetscScalar ***)((*a) + m);
3185: c = (PetscScalar **)(b + m * n);
3186: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3187: for (i = 0; i < m; i++)
3188: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3189: for (i = 0; i < m; i++)
3190: for (j = 0; j < n; j++)
3191: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3192: *a -= mstart;
3193: PetscFunctionReturn(PETSC_SUCCESS);
3194: }
3196: /*@C
3197: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3198: processor's portion of the vector data. You MUST call `VecRestoreArray4dWrite()`
3199: when you no longer need access to the array.
3201: Logically Collective
3203: Input Parameters:
3204: + x - the vector
3205: . m - first dimension of four dimensional array
3206: . n - second dimension of four dimensional array
3207: . p - third dimension of four dimensional array
3208: . q - fourth dimension of four dimensional array
3209: . mstart - first index you will use in first coordinate direction (often 0)
3210: . nstart - first index in the second coordinate direction (often 0)
3211: . pstart - first index in the third coordinate direction (often 0)
3212: - qstart - first index in the fourth coordinate direction (often 0)
3214: Output Parameter:
3215: . a - location to put pointer to the array
3217: Level: developer
3219: Notes:
3220: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3221: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3222: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3223: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3225: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3227: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3228: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3229: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3230: @*/
3231: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3232: {
3233: PetscInt i, N, j, k;
3234: PetscScalar *aa, ***b, **c;
3236: PetscFunctionBegin;
3238: PetscAssertPointer(a, 10);
3240: PetscCall(VecGetLocalSize(x, &N));
3241: 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);
3242: PetscCall(VecGetArrayWrite(x, &aa));
3244: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3245: b = (PetscScalar ***)((*a) + m);
3246: c = (PetscScalar **)(b + m * n);
3247: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3248: for (i = 0; i < m; i++)
3249: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3250: for (i = 0; i < m; i++)
3251: for (j = 0; j < n; j++)
3252: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3253: *a -= mstart;
3254: PetscFunctionReturn(PETSC_SUCCESS);
3255: }
3257: /*@C
3258: VecRestoreArray4d - Restores a vector after `VecGetArray4d()` has been called.
3260: Logically Collective
3262: Input Parameters:
3263: + x - the vector
3264: . m - first dimension of four dimensional array
3265: . n - second dimension of the four dimensional array
3266: . p - third dimension of the four dimensional array
3267: . q - fourth dimension of the four dimensional array
3268: . mstart - first index you will use in first coordinate direction (often 0)
3269: . nstart - first index in the second coordinate direction (often 0)
3270: . pstart - first index in the third coordinate direction (often 0)
3271: . qstart - first index in the fourth coordinate direction (often 0)
3272: - a - location of pointer to array obtained from VecGetArray4d()
3274: Level: developer
3276: Notes:
3277: For regular PETSc vectors this routine does not involve any copies. For
3278: any special vectors that do not store local vector data in a contiguous
3279: array, this routine will copy the data back into the underlying
3280: vector data structure from the array obtained with `VecGetArray()`.
3282: This routine actually zeros out the `a` pointer.
3284: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3285: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3286: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`
3287: @*/
3288: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3289: {
3290: void *dummy;
3292: PetscFunctionBegin;
3294: PetscAssertPointer(a, 10);
3296: dummy = (void *)(*a + mstart);
3297: PetscCall(PetscFree(dummy));
3298: PetscCall(VecRestoreArray(x, NULL));
3299: *a = NULL;
3300: PetscFunctionReturn(PETSC_SUCCESS);
3301: }
3303: /*@C
3304: VecRestoreArray4dWrite - Restores a vector after `VecGetArray4dWrite()` has been called.
3306: Logically Collective
3308: Input Parameters:
3309: + x - the vector
3310: . m - first dimension of four dimensional array
3311: . n - second dimension of the four dimensional array
3312: . p - third dimension of the four dimensional array
3313: . q - fourth dimension of the four dimensional array
3314: . mstart - first index you will use in first coordinate direction (often 0)
3315: . nstart - first index in the second coordinate direction (often 0)
3316: . pstart - first index in the third coordinate direction (often 0)
3317: . qstart - first index in the fourth coordinate direction (often 0)
3318: - a - location of pointer to array obtained from `VecGetArray4d()`
3320: Level: developer
3322: Notes:
3323: For regular PETSc vectors this routine does not involve any copies. For
3324: any special vectors that do not store local vector data in a contiguous
3325: array, this routine will copy the data back into the underlying
3326: vector data structure from the array obtained with `VecGetArray()`.
3328: This routine actually zeros out the `a` pointer.
3330: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3331: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3332: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3333: @*/
3334: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3335: {
3336: void *dummy;
3338: PetscFunctionBegin;
3340: PetscAssertPointer(a, 10);
3342: dummy = (void *)(*a + mstart);
3343: PetscCall(PetscFree(dummy));
3344: PetscCall(VecRestoreArrayWrite(x, NULL));
3345: *a = NULL;
3346: PetscFunctionReturn(PETSC_SUCCESS);
3347: }
3349: /*@C
3350: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3351: processor's portion of the vector data. You MUST call `VecRestoreArray2dRead()`
3352: when you no longer need access to the array.
3354: Logically Collective
3356: Input Parameters:
3357: + x - the vector
3358: . m - first dimension of two dimensional array
3359: . n - second dimension of two dimensional array
3360: . mstart - first index you will use in first coordinate direction (often 0)
3361: - nstart - first index in the second coordinate direction (often 0)
3363: Output Parameter:
3364: . a - location to put pointer to the array
3366: Level: developer
3368: Notes:
3369: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
3370: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3371: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3372: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
3374: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3376: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3377: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3378: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3379: @*/
3380: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3381: {
3382: PetscInt i, N;
3383: const PetscScalar *aa;
3385: PetscFunctionBegin;
3387: PetscAssertPointer(a, 6);
3389: PetscCall(VecGetLocalSize(x, &N));
3390: 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);
3391: PetscCall(VecGetArrayRead(x, &aa));
3393: PetscCall(PetscMalloc1(m, a));
3394: for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3395: *a -= mstart;
3396: PetscFunctionReturn(PETSC_SUCCESS);
3397: }
3399: /*@C
3400: VecRestoreArray2dRead - Restores a vector after `VecGetArray2dRead()` has been called.
3402: Logically Collective
3404: Input Parameters:
3405: + x - the vector
3406: . m - first dimension of two dimensional array
3407: . n - second dimension of the two dimensional array
3408: . mstart - first index you will use in first coordinate direction (often 0)
3409: . nstart - first index in the second coordinate direction (often 0)
3410: - a - location of pointer to array obtained from VecGetArray2d()
3412: Level: developer
3414: Notes:
3415: For regular PETSc vectors this routine does not involve any copies. For
3416: any special vectors that do not store local vector data in a contiguous
3417: array, this routine will copy the data back into the underlying
3418: vector data structure from the array obtained with `VecGetArray()`.
3420: This routine actually zeros out the `a` pointer.
3422: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3423: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3424: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3425: @*/
3426: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3427: {
3428: void *dummy;
3430: PetscFunctionBegin;
3432: PetscAssertPointer(a, 6);
3434: dummy = (void *)(*a + mstart);
3435: PetscCall(PetscFree(dummy));
3436: PetscCall(VecRestoreArrayRead(x, NULL));
3437: *a = NULL;
3438: PetscFunctionReturn(PETSC_SUCCESS);
3439: }
3441: /*@C
3442: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3443: processor's portion of the vector data. You MUST call `VecRestoreArray1dRead()`
3444: when you no longer need access to the array.
3446: Logically Collective
3448: Input Parameters:
3449: + x - the vector
3450: . m - first dimension of two dimensional array
3451: - mstart - first index you will use in first coordinate direction (often 0)
3453: Output Parameter:
3454: . a - location to put pointer to the array
3456: Level: developer
3458: Notes:
3459: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3460: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3461: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3463: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3465: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3466: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3467: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3468: @*/
3469: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3470: {
3471: PetscInt N;
3473: PetscFunctionBegin;
3475: PetscAssertPointer(a, 4);
3477: PetscCall(VecGetLocalSize(x, &N));
3478: 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);
3479: PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3480: *a -= mstart;
3481: PetscFunctionReturn(PETSC_SUCCESS);
3482: }
3484: /*@C
3485: VecRestoreArray1dRead - Restores a vector after `VecGetArray1dRead()` has been called.
3487: Logically Collective
3489: Input Parameters:
3490: + x - the vector
3491: . m - first dimension of two dimensional array
3492: . mstart - first index you will use in first coordinate direction (often 0)
3493: - a - location of pointer to array obtained from `VecGetArray1dRead()`
3495: Level: developer
3497: Notes:
3498: For regular PETSc vectors this routine does not involve any copies. For
3499: any special vectors that do not store local vector data in a contiguous
3500: array, this routine will copy the data back into the underlying
3501: vector data structure from the array obtained with `VecGetArray1dRead()`.
3503: This routine actually zeros out the `a` pointer.
3505: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3506: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3507: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3508: @*/
3509: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3510: {
3511: PetscFunctionBegin;
3514: PetscCall(VecRestoreArrayRead(x, NULL));
3515: *a = NULL;
3516: PetscFunctionReturn(PETSC_SUCCESS);
3517: }
3519: /*@C
3520: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3521: processor's portion of the vector data. You MUST call `VecRestoreArray3dRead()`
3522: when you no longer need access to the array.
3524: Logically Collective
3526: Input Parameters:
3527: + x - the vector
3528: . m - first dimension of three dimensional array
3529: . n - second dimension of three dimensional array
3530: . p - third dimension of three dimensional array
3531: . mstart - first index you will use in first coordinate direction (often 0)
3532: . nstart - first index in the second coordinate direction (often 0)
3533: - pstart - first index in the third coordinate direction (often 0)
3535: Output Parameter:
3536: . a - location to put pointer to the array
3538: Level: developer
3540: Notes:
3541: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3542: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3543: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3544: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3dRead()`.
3546: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3548: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3549: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3550: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3551: @*/
3552: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3553: {
3554: PetscInt i, N, j;
3555: const PetscScalar *aa;
3556: PetscScalar **b;
3558: PetscFunctionBegin;
3560: PetscAssertPointer(a, 8);
3562: PetscCall(VecGetLocalSize(x, &N));
3563: 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);
3564: PetscCall(VecGetArrayRead(x, &aa));
3566: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3567: b = (PetscScalar **)((*a) + m);
3568: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3569: for (i = 0; i < m; i++)
3570: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset((PetscScalar *)aa, i * n * p + j * p - pstart);
3571: *a -= mstart;
3572: PetscFunctionReturn(PETSC_SUCCESS);
3573: }
3575: /*@C
3576: VecRestoreArray3dRead - Restores a vector after `VecGetArray3dRead()` has been called.
3578: Logically Collective
3580: Input Parameters:
3581: + x - the vector
3582: . m - first dimension of three dimensional array
3583: . n - second dimension of the three dimensional array
3584: . p - third dimension of the three dimensional array
3585: . mstart - first index you will use in first coordinate direction (often 0)
3586: . nstart - first index in the second coordinate direction (often 0)
3587: . pstart - first index in the third coordinate direction (often 0)
3588: - a - location of pointer to array obtained from `VecGetArray3dRead()`
3590: Level: developer
3592: Notes:
3593: For regular PETSc vectors this routine does not involve any copies. For
3594: any special vectors that do not store local vector data in a contiguous
3595: array, this routine will copy the data back into the underlying
3596: vector data structure from the array obtained with `VecGetArray()`.
3598: This routine actually zeros out the `a` pointer.
3600: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3601: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3602: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3603: @*/
3604: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3605: {
3606: void *dummy;
3608: PetscFunctionBegin;
3610: PetscAssertPointer(a, 8);
3612: dummy = (void *)(*a + mstart);
3613: PetscCall(PetscFree(dummy));
3614: PetscCall(VecRestoreArrayRead(x, NULL));
3615: *a = NULL;
3616: PetscFunctionReturn(PETSC_SUCCESS);
3617: }
3619: /*@C
3620: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3621: processor's portion of the vector data. You MUST call `VecRestoreArray4dRead()`
3622: when you no longer need access to the array.
3624: Logically Collective
3626: Input Parameters:
3627: + x - the vector
3628: . m - first dimension of four dimensional array
3629: . n - second dimension of four dimensional array
3630: . p - third dimension of four dimensional array
3631: . q - fourth dimension of four dimensional array
3632: . mstart - first index you will use in first coordinate direction (often 0)
3633: . nstart - first index in the second coordinate direction (often 0)
3634: . pstart - first index in the third coordinate direction (often 0)
3635: - qstart - first index in the fourth coordinate direction (often 0)
3637: Output Parameter:
3638: . a - location to put pointer to the array
3640: Level: beginner
3642: Notes:
3643: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3644: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3645: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3646: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3648: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3650: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3651: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3652: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3653: @*/
3654: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3655: {
3656: PetscInt i, N, j, k;
3657: const PetscScalar *aa;
3658: PetscScalar ***b, **c;
3660: PetscFunctionBegin;
3662: PetscAssertPointer(a, 10);
3664: PetscCall(VecGetLocalSize(x, &N));
3665: 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);
3666: PetscCall(VecGetArrayRead(x, &aa));
3668: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3669: b = (PetscScalar ***)((*a) + m);
3670: c = (PetscScalar **)(b + m * n);
3671: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3672: for (i = 0; i < m; i++)
3673: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3674: for (i = 0; i < m; i++)
3675: for (j = 0; j < n; j++)
3676: 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;
3677: *a -= mstart;
3678: PetscFunctionReturn(PETSC_SUCCESS);
3679: }
3681: /*@C
3682: VecRestoreArray4dRead - Restores a vector after `VecGetArray4d()` has been called.
3684: Logically Collective
3686: Input Parameters:
3687: + x - the vector
3688: . m - first dimension of four dimensional array
3689: . n - second dimension of the four dimensional array
3690: . p - third dimension of the four dimensional array
3691: . q - fourth dimension of the four dimensional array
3692: . mstart - first index you will use in first coordinate direction (often 0)
3693: . nstart - first index in the second coordinate direction (often 0)
3694: . pstart - first index in the third coordinate direction (often 0)
3695: . qstart - first index in the fourth coordinate direction (often 0)
3696: - a - location of pointer to array obtained from `VecGetArray4dRead()`
3698: Level: beginner
3700: Notes:
3701: For regular PETSc vectors this routine does not involve any copies. For
3702: any special vectors that do not store local vector data in a contiguous
3703: array, this routine will copy the data back into the underlying
3704: vector data structure from the array obtained with `VecGetArray()`.
3706: This routine actually zeros out the `a` pointer.
3708: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3709: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3710: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3711: @*/
3712: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3713: {
3714: void *dummy;
3716: PetscFunctionBegin;
3718: PetscAssertPointer(a, 10);
3720: dummy = (void *)(*a + mstart);
3721: PetscCall(PetscFree(dummy));
3722: PetscCall(VecRestoreArrayRead(x, NULL));
3723: *a = NULL;
3724: PetscFunctionReturn(PETSC_SUCCESS);
3725: }
3727: /*@
3728: VecLockGet - Get the current lock status of a vector
3730: Logically Collective
3732: Input Parameter:
3733: . x - the vector
3735: Output Parameter:
3736: . state - greater than zero indicates the vector is locked for read; less than zero indicates the vector is
3737: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3739: Level: advanced
3741: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3742: @*/
3743: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3744: {
3745: PetscFunctionBegin;
3747: PetscAssertPointer(state, 2);
3748: *state = x->lock;
3749: PetscFunctionReturn(PETSC_SUCCESS);
3750: }
3752: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3753: {
3754: PetscFunctionBegin;
3756: PetscAssertPointer(file, 2);
3757: PetscAssertPointer(func, 3);
3758: PetscAssertPointer(line, 4);
3759: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3760: {
3761: const int index = x->lockstack.currentsize - 1;
3763: *file = index < 0 ? NULL : x->lockstack.file[index];
3764: *func = index < 0 ? NULL : x->lockstack.function[index];
3765: *line = index < 0 ? 0 : x->lockstack.line[index];
3766: }
3767: #else
3768: *file = NULL;
3769: *func = NULL;
3770: *line = 0;
3771: #endif
3772: PetscFunctionReturn(PETSC_SUCCESS);
3773: }
3775: /*@
3776: VecLockReadPush - Push a read-only lock on a vector to prevent it from being written to
3778: Logically Collective
3780: Input Parameter:
3781: . x - the vector
3783: Level: intermediate
3785: Notes:
3786: If this is set then calls to `VecGetArray()` or `VecSetValues()` or any other routines that change the vectors values will generate an error.
3788: The call can be nested, i.e., called multiple times on the same vector, but each `VecLockReadPush()` has to have one matching
3789: `VecLockReadPop()`, which removes the latest read-only lock.
3791: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3792: @*/
3793: PetscErrorCode VecLockReadPush(Vec x)
3794: {
3795: PetscFunctionBegin;
3797: 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");
3798: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3799: {
3800: const char *file, *func;
3801: int index, line;
3803: if ((index = petscstack.currentsize - 2) < 0) {
3804: // vec was locked "outside" of petsc, either in user-land or main. the error message will
3805: // now show this function as the culprit, but it will include the stacktrace
3806: file = "unknown user-file";
3807: func = "unknown_user_function";
3808: line = 0;
3809: } else {
3810: file = petscstack.file[index];
3811: func = petscstack.function[index];
3812: line = petscstack.line[index];
3813: }
3814: PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
3815: }
3816: #endif
3817: PetscFunctionReturn(PETSC_SUCCESS);
3818: }
3820: /*@
3821: VecLockReadPop - Pop a read-only lock from a vector
3823: Logically Collective
3825: Input Parameter:
3826: . x - the vector
3828: Level: intermediate
3830: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
3831: @*/
3832: PetscErrorCode VecLockReadPop(Vec x)
3833: {
3834: PetscFunctionBegin;
3836: PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
3837: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3838: {
3839: const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];
3841: PetscStackPop_Private(x->lockstack, previous);
3842: }
3843: #endif
3844: PetscFunctionReturn(PETSC_SUCCESS);
3845: }
3847: /*@
3848: VecLockWriteSet - Lock or unlock a vector for exclusive read/write access
3850: Logically Collective
3852: Input Parameters:
3853: + x - the vector
3854: - flg - `PETSC_TRUE` to lock the vector for exclusive read/write access; `PETSC_FALSE` to unlock it.
3856: Level: intermediate
3858: Notes:
3859: The function is useful in split-phase computations, which usually have a begin phase and an end phase.
3860: One can call `VecLockWriteSet`(x,`PETSC_TRUE`) in the begin phase to lock a vector for exclusive
3861: access, and call `VecLockWriteSet`(x,`PETSC_FALSE`) in the end phase to unlock the vector from exclusive
3862: access. In this way, one is ensured no other operations can access the vector in between. The code may like
3864: .vb
3865: VecGetArray(x,&xdata); // begin phase
3866: VecLockWriteSet(v,PETSC_TRUE);
3868: Other operations, which can not access x anymore (they can access xdata, of course)
3870: VecRestoreArray(x,&vdata); // end phase
3871: VecLockWriteSet(v,PETSC_FALSE);
3872: .ve
3874: The call can not be nested on the same vector, in other words, one can not call `VecLockWriteSet`(x,`PETSC_TRUE`)
3875: again before calling `VecLockWriteSet`(v,`PETSC_FALSE`).
3877: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
3878: @*/
3879: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
3880: {
3881: PetscFunctionBegin;
3883: if (flg) {
3884: 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");
3885: 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");
3886: x->lock = -1;
3887: } else {
3888: 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");
3889: x->lock = 0;
3890: }
3891: PetscFunctionReturn(PETSC_SUCCESS);
3892: }