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: Note:
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;
204: PetscCall(VecLockReadPush(x));
207: VecCheckAssembled(x);
209: PetscAssertPointer(val, 3);
211: PetscCall(VecNormAvailable(x, type, &flg, val));
212: // check that all MPI processes call this routine together and have same availability
213: if (PetscDefined(USE_DEBUG)) {
214: PetscMPIInt b0 = (PetscMPIInt)flg, b2[2];
215: b2[0] = -b0;
216: b2[1] = b0;
217: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, b2, 2, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
218: 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]);
219: if (flg) {
220: PetscReal b2[2];
221: b2[0] = -(*val);
222: b2[1] = *val;
223: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, b2, 2, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)x)));
224: 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);
225: }
226: }
227: if (!flg) {
228: PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
229: PetscUseTypeMethod(x, norm, type, val);
230: PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
232: if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
233: }
234: PetscCall(VecLockReadPop(x));
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: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
256: `VecNormBegin()`, `VecNormEnd()`
257: @*/
258: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
259: {
260: PetscFunctionBegin;
263: PetscAssertPointer(available, 3);
264: PetscAssertPointer(val, 4);
266: if (type == NORM_1_AND_2) {
267: *available = PETSC_FALSE;
268: } else {
269: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
270: }
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: /*@
275: VecNormalize - Normalizes a vector by its 2-norm.
277: Collective
279: Input Parameter:
280: . x - the vector
282: Output Parameter:
283: . val - the vector norm before normalization. May be `NULL` if the value is not needed.
285: Level: intermediate
287: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
288: @*/
289: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
290: {
291: PetscReal norm;
293: PetscFunctionBegin;
296: PetscCall(VecSetErrorIfLocked(x, 1));
297: if (val) PetscAssertPointer(val, 2);
298: PetscCall(PetscLogEventBegin(VEC_Normalize, x, 0, 0, 0));
299: PetscCall(VecNorm(x, NORM_2, &norm));
300: if (norm == 0.0) PetscCall(PetscInfo(x, "Vector of zero norm can not be normalized; Returning only the zero norm\n"));
301: else if (PetscIsInfOrNanReal(norm)) PetscCall(PetscInfo(x, "Vector with infinity or NaN norm can not be normalized; Returning only the norm\n"));
302: else {
303: PetscScalar s = 1.0 / norm;
304: PetscCall(VecScale(x, s));
305: }
306: PetscCall(PetscLogEventEnd(VEC_Normalize, x, 0, 0, 0));
307: if (val) *val = norm;
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: /*@
312: VecMax - Determines the vector component with maximum real part and its location.
314: Collective
316: Input Parameter:
317: . x - the vector
319: Output Parameters:
320: + p - the index of `val` (pass `NULL` if you don't want this) in the vector
321: - val - the maximum component
323: Level: intermediate
325: Notes:
326: Returns the value `PETSC_MIN_REAL` and negative `p` if the vector is of length 0.
328: Returns the smallest index with the maximum value
330: Developer Note:
331: The Nag Fortran compiler does not like the symbol name VecMax
333: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
334: @*/
335: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
336: {
337: PetscFunctionBegin;
340: VecCheckAssembled(x);
341: if (p) PetscAssertPointer(p, 2);
342: PetscAssertPointer(val, 3);
343: PetscCall(VecLockReadPush(x));
344: PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
345: PetscUseTypeMethod(x, max, p, val);
346: PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
347: PetscCall(VecLockReadPop(x));
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: /*@
352: VecMin - Determines the vector component with minimum real part and its location.
354: Collective
356: Input Parameter:
357: . x - the vector
359: Output Parameters:
360: + p - the index of `val` (pass `NULL` if you don't want this location) in the vector
361: - val - the minimum component
363: Level: intermediate
365: Notes:
366: Returns the value `PETSC_MAX_REAL` and negative `p` if the vector is of length 0.
368: This returns the smallest index with the minimum value
370: Developer Note:
371: The Nag Fortran compiler does not like the symbol name VecMin
373: .seealso: [](ch_vectors), `Vec`, `VecMax()`
374: @*/
375: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
376: {
377: PetscFunctionBegin;
380: VecCheckAssembled(x);
381: if (p) PetscAssertPointer(p, 2);
382: PetscAssertPointer(val, 3);
383: PetscCall(VecLockReadPush(x));
384: PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
385: PetscUseTypeMethod(x, min, p, val);
386: PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
387: PetscCall(VecLockReadPop(x));
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: /*@
392: VecTDot - Computes an indefinite vector dot product. That is, this
393: routine does NOT use the complex conjugate.
395: Collective
397: Input Parameters:
398: + x - first vector
399: - y - second vector
401: Output Parameter:
402: . val - the dot product
404: Level: intermediate
406: Notes for Users of Complex Numbers:
407: For complex vectors, `VecTDot()` computes the indefinite form
408: .vb
409: val = (x,y) = y^T x,
410: .ve
411: where y^T denotes the transpose of y.
413: Use `VecDot()` for the inner product
414: .vb
415: val = (x,y) = y^H x,
416: .ve
417: where y^H denotes the conjugate transpose of y.
419: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
420: @*/
421: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
422: {
423: PetscFunctionBegin;
426: PetscAssertPointer(val, 3);
429: PetscCheckSameTypeAndComm(x, 1, y, 2);
430: VecCheckSameSize(x, 1, y, 2);
431: VecCheckAssembled(x);
432: VecCheckAssembled(y);
434: PetscCall(VecLockReadPush(x));
435: PetscCall(VecLockReadPush(y));
436: PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
437: PetscUseTypeMethod(x, tdot, y, val);
438: PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
439: PetscCall(VecLockReadPop(x));
440: PetscCall(VecLockReadPop(y));
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
444: PetscErrorCode VecScaleAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
445: {
446: PetscReal norms[4];
447: PetscBool flgs[4];
448: PetscScalar one = 1.0;
450: PetscFunctionBegin;
453: VecCheckAssembled(x);
454: PetscCall(VecSetErrorIfLocked(x, 1));
456: if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);
458: /* get current stashed norms */
459: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
461: PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
462: VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
463: PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));
465: PetscCall(PetscObjectStateIncrease((PetscObject)x));
466: /* put the scaled stashed norms back into the Vec */
467: for (PetscInt i = 0; i < 4; i++) {
468: PetscReal ar = PetscAbsScalar(alpha);
469: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
470: }
471: PetscFunctionReturn(PETSC_SUCCESS);
472: }
474: /*@
475: VecScale - Scales a vector.
477: Logically Collective
479: Input Parameters:
480: + x - the vector
481: - alpha - the scalar
483: Level: intermediate
485: Note:
486: For a vector with n components, `VecScale()` computes x[i] = alpha * x[i], for i=1,...,n.
488: .seealso: [](ch_vectors), `Vec`, `VecSet()`
489: @*/
490: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
491: {
492: PetscFunctionBegin;
493: PetscCall(VecScaleAsync_Private(x, alpha, NULL));
494: PetscFunctionReturn(PETSC_SUCCESS);
495: }
497: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
498: {
499: PetscFunctionBegin;
502: VecCheckAssembled(x);
504: PetscCall(VecSetErrorIfLocked(x, 1));
506: if (alpha == 0) {
507: PetscReal norm;
508: PetscBool set;
510: PetscCall(VecNormAvailable(x, NORM_2, &set, &norm));
511: if (set == PETSC_TRUE && norm == 0) PetscFunctionReturn(PETSC_SUCCESS);
512: }
513: PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
514: VecMethodDispatch(x, dctx, VecAsyncFnName(Set), set, (Vec, PetscScalar, PetscDeviceContext), alpha);
515: PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
516: PetscCall(PetscObjectStateIncrease((PetscObject)x));
518: /* norms can be simply set (if |alpha|*N not too large) */
519: {
520: PetscReal val = PetscAbsScalar(alpha);
521: const PetscInt N = x->map->N;
523: if (N == 0) {
524: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0));
525: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
526: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
527: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
528: } else if (val > PETSC_MAX_REAL / N) {
529: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
530: } else {
531: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
532: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
533: val *= PetscSqrtReal((PetscReal)N);
534: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
535: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
536: }
537: }
538: PetscFunctionReturn(PETSC_SUCCESS);
539: }
541: /*@
542: VecSet - Sets all components of a vector to a single scalar value.
544: Logically Collective
546: Input Parameters:
547: + x - the vector
548: - alpha - the scalar
550: Level: beginner
552: Notes:
553: For a vector of dimension n, `VecSet()` sets x[i] = alpha, for i=1,...,n,
554: so that all vector entries then equal the identical
555: scalar value, `alpha`. Use the more general routine
556: `VecSetValues()` to set different vector entries.
558: You CANNOT call this after you have called `VecSetValues()` but before you call
559: `VecAssemblyBegin()`
561: If `alpha` is zero and the norm of the vector is known to be zero then this skips the unneeded zeroing process
563: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
564: @*/
565: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
566: {
567: PetscFunctionBegin;
568: PetscCall(VecSetAsync_Private(x, alpha, NULL));
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: /*@
573: VecSetStdBasis - Set the vector to the i-th standard basis vector
575: Logically Collective
577: Input Parameters:
578: + x - the vector
579: - i - the component to be set to one
581: Level: beginner
583: Note:
584: This function sets x[i] = 1, and 0 otherwise.
586: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`
587: @*/
588: PetscErrorCode VecSetStdBasis(Vec x, PetscInt i)
589: {
590: PetscFunctionBegin;
593: PetscCall(VecSetErrorIfLocked(x, 1));
594: if (x->ops->setstdbasis) PetscUseTypeMethod(x, setstdbasis, i);
595: else {
596: PetscInt st, en;
598: PetscCall(VecGetOwnershipRange(x, &st, &en));
599: PetscCall(VecSet(x, 0.));
600: if (st <= i && i < en) PetscCall(VecSetValue(x, i, 1.0, INSERT_VALUES));
601: PetscCall(VecAssemblyBegin(x));
602: PetscCall(VecAssemblyEnd(x));
603: }
604: PetscCall(PetscObjectStateIncrease((PetscObject)x));
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }
608: PetscErrorCode VecAXPYAsync_Private(Vec y, PetscScalar alpha, Vec x, PetscDeviceContext dctx)
609: {
610: PetscFunctionBegin;
615: PetscCheckSameTypeAndComm(x, 3, y, 1);
616: VecCheckSameSize(x, 3, y, 1);
617: VecCheckAssembled(x);
618: VecCheckAssembled(y);
620: if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
621: PetscCall(VecSetErrorIfLocked(y, 1));
622: if (x == y) {
623: PetscCall(VecScale(y, alpha + 1.0));
624: PetscFunctionReturn(PETSC_SUCCESS);
625: }
626: PetscCall(VecLockReadPush(x));
627: PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
628: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPY), axpy, (Vec, PetscScalar, Vec, PetscDeviceContext), alpha, x);
629: PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
630: PetscCall(VecLockReadPop(x));
631: PetscCall(PetscObjectStateIncrease((PetscObject)y));
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }
635: /*@
636: VecAXPY - Computes `y = alpha x + y`.
638: Logically Collective
640: Input Parameters:
641: + alpha - the scalar
642: . x - vector scale by `alpha`
643: - y - vector accumulated into
645: Output Parameter:
646: . y - output vector
648: Level: intermediate
650: Notes:
651: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
652: .vb
653: VecAXPY(y,alpha,x) y = alpha x + y
654: VecAYPX(y,beta,x) y = x + beta y
655: VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
656: VecWAXPY(w,alpha,x,y) w = alpha x + y
657: VecAXPBYPCZ(z,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
658: VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
659: .ve
661: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
662: @*/
663: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
664: {
665: PetscFunctionBegin;
666: PetscCall(VecAXPYAsync_Private(y, alpha, x, NULL));
667: PetscFunctionReturn(PETSC_SUCCESS);
668: }
670: PetscErrorCode VecAYPXAsync_Private(Vec y, PetscScalar beta, Vec x, PetscDeviceContext dctx)
671: {
672: PetscFunctionBegin;
677: PetscCheckSameTypeAndComm(x, 3, y, 1);
678: VecCheckSameSize(x, 1, y, 3);
679: VecCheckAssembled(x);
680: VecCheckAssembled(y);
682: PetscCall(VecSetErrorIfLocked(y, 1));
683: if (x == y) {
684: PetscCall(VecScale(y, beta + 1.0));
685: PetscFunctionReturn(PETSC_SUCCESS);
686: }
687: PetscCall(VecLockReadPush(x));
688: if (beta == (PetscScalar)0.0) {
689: PetscCall(VecCopy(x, y));
690: } else {
691: PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
692: VecMethodDispatch(y, dctx, VecAsyncFnName(AYPX), aypx, (Vec, PetscScalar, Vec, PetscDeviceContext), beta, x);
693: PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
694: PetscCall(PetscObjectStateIncrease((PetscObject)y));
695: }
696: PetscCall(VecLockReadPop(x));
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: /*@
701: VecAYPX - Computes `y = x + beta y`.
703: Logically Collective
705: Input Parameters:
706: + beta - the scalar
707: . x - the unscaled vector
708: - y - the vector to be scaled
710: Output Parameter:
711: . y - output vector
713: Level: intermediate
715: Developer Notes:
716: The implementation is optimized for `beta` of -1.0, 0.0, and 1.0
718: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
719: @*/
720: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
721: {
722: PetscFunctionBegin;
723: PetscCall(VecAYPXAsync_Private(y, beta, x, NULL));
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: PetscErrorCode VecAXPBYAsync_Private(Vec y, PetscScalar alpha, PetscScalar beta, Vec x, PetscDeviceContext dctx)
728: {
729: PetscFunctionBegin;
734: PetscCheckSameTypeAndComm(x, 4, y, 1);
735: VecCheckSameSize(y, 1, x, 4);
736: VecCheckAssembled(x);
737: VecCheckAssembled(y);
740: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
741: if (x == y) {
742: PetscCall(VecScale(y, alpha + beta));
743: PetscFunctionReturn(PETSC_SUCCESS);
744: }
746: PetscCall(VecSetErrorIfLocked(y, 1));
747: PetscCall(VecLockReadPush(x));
748: PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
749: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPBY), axpby, (Vec, PetscScalar, PetscScalar, Vec, PetscDeviceContext), alpha, beta, x);
750: PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
751: PetscCall(PetscObjectStateIncrease((PetscObject)y));
752: PetscCall(VecLockReadPop(x));
753: PetscFunctionReturn(PETSC_SUCCESS);
754: }
756: /*@
757: VecAXPBY - Computes `y = alpha x + beta y`.
759: Logically Collective
761: Input Parameters:
762: + alpha - first scalar
763: . beta - second scalar
764: . x - the first scaled vector
765: - y - the second scaled vector
767: Output Parameter:
768: . y - output vector
770: Level: intermediate
772: Developer Notes:
773: The implementation is optimized for `alpha` and/or `beta` values of 0.0 and 1.0
775: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
776: @*/
777: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
778: {
779: PetscFunctionBegin;
780: PetscCall(VecAXPBYAsync_Private(y, alpha, beta, x, NULL));
781: PetscFunctionReturn(PETSC_SUCCESS);
782: }
784: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
785: {
786: PetscFunctionBegin;
793: PetscCheckSameTypeAndComm(x, 5, y, 6);
794: PetscCheckSameTypeAndComm(x, 5, z, 1);
795: VecCheckSameSize(x, 5, y, 6);
796: VecCheckSameSize(x, 5, z, 1);
797: PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
798: PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
799: VecCheckAssembled(x);
800: VecCheckAssembled(y);
801: VecCheckAssembled(z);
805: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
807: PetscCall(VecSetErrorIfLocked(z, 1));
808: PetscCall(VecLockReadPush(x));
809: PetscCall(VecLockReadPush(y));
810: PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
811: VecMethodDispatch(z, dctx, VecAsyncFnName(AXPBYPCZ), axpbypcz, (Vec, PetscScalar, PetscScalar, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, beta, gamma, x, y);
812: PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
813: PetscCall(PetscObjectStateIncrease((PetscObject)z));
814: PetscCall(VecLockReadPop(x));
815: PetscCall(VecLockReadPop(y));
816: PetscFunctionReturn(PETSC_SUCCESS);
817: }
818: /*@
819: VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`
821: Logically Collective
823: Input Parameters:
824: + alpha - first scalar
825: . beta - second scalar
826: . gamma - third scalar
827: . x - first vector
828: . y - second vector
829: - z - third vector
831: Output Parameter:
832: . z - output vector
834: Level: intermediate
836: Note:
837: `x`, `y` and `z` must be different vectors
839: Developer Notes:
840: The implementation is optimized for `alpha` of 1.0 and `gamma` of 1.0 or 0.0
842: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
843: @*/
844: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
845: {
846: PetscFunctionBegin;
847: PetscCall(VecAXPBYPCZAsync_Private(z, alpha, beta, gamma, x, y, NULL));
848: PetscFunctionReturn(PETSC_SUCCESS);
849: }
851: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
852: {
853: PetscFunctionBegin;
860: PetscCheckSameTypeAndComm(x, 3, y, 4);
861: PetscCheckSameTypeAndComm(y, 4, w, 1);
862: VecCheckSameSize(x, 3, y, 4);
863: VecCheckSameSize(x, 3, w, 1);
864: PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
865: PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
866: VecCheckAssembled(x);
867: VecCheckAssembled(y);
869: PetscCall(VecSetErrorIfLocked(w, 1));
871: PetscCall(VecLockReadPush(x));
872: PetscCall(VecLockReadPush(y));
873: if (alpha == (PetscScalar)0.0) {
874: PetscCall(VecCopyAsync_Private(y, w, dctx));
875: } else {
876: PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
877: VecMethodDispatch(w, dctx, VecAsyncFnName(WAXPY), waxpy, (Vec, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, x, y);
878: PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
879: PetscCall(PetscObjectStateIncrease((PetscObject)w));
880: }
881: PetscCall(VecLockReadPop(x));
882: PetscCall(VecLockReadPop(y));
883: PetscFunctionReturn(PETSC_SUCCESS);
884: }
886: /*@
887: VecWAXPY - Computes `w = alpha x + y`.
889: Logically Collective
891: Input Parameters:
892: + alpha - the scalar
893: . x - first vector, multiplied by `alpha`
894: - y - second vector
896: Output Parameter:
897: . w - the result
899: Level: intermediate
901: Note:
902: `w` cannot be either `x` or `y`, but `x` and `y` can be the same
904: Developer Notes:
905: The implementation is optimized for alpha of -1.0, 0.0, and 1.0
907: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
908: @*/
909: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
910: {
911: PetscFunctionBegin;
912: PetscCall(VecWAXPYAsync_Private(w, alpha, x, y, NULL));
913: PetscFunctionReturn(PETSC_SUCCESS);
914: }
916: /*@
917: VecSetValues - Inserts or adds values into certain locations of a vector.
919: Not Collective
921: Input Parameters:
922: + x - vector to insert in
923: . ni - number of elements to add
924: . ix - indices where to add
925: . y - array of values. Pass `NULL` to set all zeroes.
926: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries
928: Level: beginner
930: Notes:
931: .vb
932: `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
933: .ve
935: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
936: options cannot be mixed without intervening calls to the assembly
937: routines.
939: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
940: MUST be called after all calls to `VecSetValues()` have been completed.
942: VecSetValues() uses 0-based indices in Fortran as well as in C.
944: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
945: negative indices may be passed in ix. These rows are
946: simply ignored. This allows easily inserting element load matrices
947: with homogeneous Dirichlet boundary conditions that you don't want represented
948: in the vector.
950: Fortran Note:
951: If any of `ix` and `y` are scalars pass them using, for example,
952: .vb
953: call VecSetValues(mat, one, [ix], [y], INSERT_VALUES, ierr)
954: .ve
956: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
957: `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`,
958: `VecOption`, `VecSetOption()`
959: @*/
960: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
961: {
962: PetscFunctionBeginHot;
964: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
965: PetscAssertPointer(ix, 3);
966: if (y) PetscAssertPointer(y, 4);
969: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
970: PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
971: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
972: PetscCall(PetscObjectStateIncrease((PetscObject)x));
973: PetscFunctionReturn(PETSC_SUCCESS);
974: }
976: /*@
977: VecGetValues - Gets values from certain locations of a vector. Currently
978: can only get values on the same processor on which they are owned
980: Not Collective
982: Input Parameters:
983: + x - vector to get values from
984: . ni - number of elements to get
985: - ix - indices where to get them from (in global 1d numbering)
987: Output Parameter:
988: . y - array of values, must be passed in with a length of `ni`
990: Level: beginner
992: Notes:
993: The user provides the allocated array y; it is NOT allocated in this routine
995: `VecGetValues()` gets y[i] = x[ix[i]], for i=0,...,ni-1.
997: `VecAssemblyBegin()` and `VecAssemblyEnd()` MUST be called before calling this if `VecSetValues()` or related routine has been called
999: VecGetValues() uses 0-based indices in Fortran as well as in C.
1001: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
1002: negative indices may be passed in ix. These rows are
1003: simply ignored.
1005: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
1006: @*/
1007: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
1008: {
1009: PetscFunctionBegin;
1011: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1012: PetscAssertPointer(ix, 3);
1013: PetscAssertPointer(y, 4);
1015: VecCheckAssembled(x);
1016: PetscUseTypeMethod(x, getvalues, ni, ix, y);
1017: PetscFunctionReturn(PETSC_SUCCESS);
1018: }
1020: /*@
1021: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
1023: Not Collective
1025: Input Parameters:
1026: + x - vector to insert in
1027: . ni - number of blocks to add
1028: . ix - indices where to add in block count, rather than element count
1029: . y - array of values. Pass `NULL` to set all zeroes.
1030: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries
1032: Level: intermediate
1034: Notes:
1035: `VecSetValuesBlocked()` sets x[bs*ix[i]+j] = y[bs*i+j],
1036: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
1038: Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
1039: options cannot be mixed without intervening calls to the assembly
1040: routines.
1042: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1043: MUST be called after all calls to `VecSetValuesBlocked()` have been completed.
1045: `VecSetValuesBlocked()` uses 0-based indices in Fortran as well as in C.
1047: Negative indices may be passed in ix, these rows are
1048: simply ignored. This allows easily inserting element load matrices
1049: with homogeneous Dirichlet boundary conditions that you don't want represented
1050: in the vector.
1052: Fortran Note:
1053: If any of `ix` and `y` are scalars pass them using, for example,
1054: .vb
1055: call VecSetValuesBlocked(mat, one, [ix], [y], INSERT_VALUES, ierr)
1056: .ve
1058: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
1059: `VecSetValues()`
1060: @*/
1061: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1062: {
1063: PetscFunctionBeginHot;
1065: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1066: PetscAssertPointer(ix, 3);
1067: if (y) PetscAssertPointer(y, 4);
1070: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1071: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1072: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1073: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1074: PetscFunctionReturn(PETSC_SUCCESS);
1075: }
1077: /*@
1078: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1079: using a local ordering of the nodes.
1081: Not Collective
1083: Input Parameters:
1084: + x - vector to insert in
1085: . ni - number of elements to add
1086: . ix - indices where to add
1087: . y - array of values. Pass `NULL` to set all zeroes.
1088: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1090: Level: intermediate
1092: Notes:
1093: `VecSetValuesLocal()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
1095: Calls to `VecSetValuesLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1096: options cannot be mixed without intervening calls to the assembly
1097: routines.
1099: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1100: MUST be called after all calls to `VecSetValuesLocal()` have been completed.
1102: `VecSetValuesLocal()` uses 0-based indices in Fortran as well as in C.
1104: Fortran Note:
1105: If any of `ix` and `y` are scalars pass them using, for example,
1106: .vb
1107: call VecSetValuesLocal(mat, one, [ix], [y], INSERT_VALUES, ierr)
1108: .ve
1110: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1111: `VecSetValuesBlockedLocal()`
1112: @*/
1113: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1114: {
1115: PetscInt lixp[128], *lix = lixp;
1117: PetscFunctionBeginHot;
1119: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1120: PetscAssertPointer(ix, 3);
1121: if (y) PetscAssertPointer(y, 4);
1124: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1125: if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
1126: if (x->map->mapping) {
1127: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1128: PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1129: PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1130: if (ni > 128) PetscCall(PetscFree(lix));
1131: } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1132: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1133: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1134: PetscFunctionReturn(PETSC_SUCCESS);
1135: }
1137: /*@
1138: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1139: using a local ordering of the nodes.
1141: Not Collective
1143: Input Parameters:
1144: + x - vector to insert in
1145: . ni - number of blocks to add
1146: . ix - indices where to add in block count, not element count
1147: . y - array of values. Pass `NULL` to set all zeroes.
1148: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1150: Level: intermediate
1152: Notes:
1153: `VecSetValuesBlockedLocal()` sets x[bs*ix[i]+j] = y[bs*i+j],
1154: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with `VecSetBlockSize()`.
1156: Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1157: options cannot be mixed without intervening calls to the assembly
1158: routines.
1160: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1161: MUST be called after all calls to `VecSetValuesBlockedLocal()` have been completed.
1163: `VecSetValuesBlockedLocal()` uses 0-based indices in Fortran as well as in C.
1165: Fortran Note:
1166: If any of `ix` and `y` are scalars pass them using, for example,
1167: .vb
1168: call VecSetValuesBlockedLocal(mat, one, [ix], [y], INSERT_VALUES, ierr)
1169: .ve
1171: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1172: `VecSetLocalToGlobalMapping()`
1173: @*/
1174: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1175: {
1176: PetscInt lixp[128], *lix = lixp;
1178: PetscFunctionBeginHot;
1180: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1181: PetscAssertPointer(ix, 3);
1182: if (y) PetscAssertPointer(y, 4);
1184: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1185: if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
1186: if (x->map->mapping) {
1187: if (ni > (PetscInt)PETSC_STATIC_ARRAY_LENGTH(lixp)) PetscCall(PetscMalloc1(ni, &lix));
1188: PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1189: PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1190: if (ni > (PetscInt)PETSC_STATIC_ARRAY_LENGTH(lixp)) PetscCall(PetscFree(lix));
1191: } else {
1192: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1193: }
1194: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1195: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1200: {
1201: PetscFunctionBegin;
1204: VecCheckAssembled(x);
1206: if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1207: PetscAssertPointer(y, 3);
1208: for (PetscInt i = 0; i < nv; ++i) {
1211: PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1212: VecCheckSameSize(x, 1, y[i], 3);
1213: VecCheckAssembled(y[i]);
1214: PetscCall(VecLockReadPush(y[i]));
1215: }
1216: PetscAssertPointer(result, 4);
1219: PetscCall(VecLockReadPush(x));
1220: PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1221: PetscCall((*mxdot)(x, nv, y, result));
1222: PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1223: PetscCall(VecLockReadPop(x));
1224: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1225: PetscFunctionReturn(PETSC_SUCCESS);
1226: }
1228: /*@
1229: VecMTDot - Computes indefinite vector multiple dot products.
1230: That is, it does NOT use the complex conjugate.
1232: Collective
1234: Input Parameters:
1235: + x - one vector
1236: . nv - number of vectors
1237: - y - array of vectors. Note that vectors are pointers
1239: Output Parameter:
1240: . val - array of the dot products
1242: Level: intermediate
1244: Notes for Users of Complex Numbers:
1245: For complex vectors, `VecMTDot()` computes the indefinite form
1246: .vb
1247: val = (x,y) = y^T x,
1248: .ve
1249: where y^T denotes the transpose of y.
1251: Use `VecMDot()` for the inner product
1252: .vb
1253: val = (x,y) = y^H x,
1254: .ve
1255: where y^H denotes the conjugate transpose of y.
1257: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1258: @*/
1259: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1260: {
1261: PetscFunctionBegin;
1263: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1264: PetscFunctionReturn(PETSC_SUCCESS);
1265: }
1267: /*@
1268: VecMDot - Computes multiple vector dot products.
1270: Collective
1272: Input Parameters:
1273: + x - one vector
1274: . nv - number of vectors
1275: - y - array of vectors.
1277: Output Parameter:
1278: . val - array of the dot products (does not allocate the array)
1280: Level: intermediate
1282: Notes for Users of Complex Numbers:
1283: For complex vectors, `VecMDot()` computes
1284: .vb
1285: val = (x,y) = y^H x,
1286: .ve
1287: where y^H denotes the conjugate transpose of y.
1289: Use `VecMTDot()` for the indefinite form
1290: .vb
1291: val = (x,y) = y^T x,
1292: .ve
1293: where y^T denotes the transpose of y.
1295: Note:
1296: The implementation may use BLAS 2 operations when the vectors `y` have been obtained with `VecDuplicateVecs()`
1298: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`, `VecDuplicateVecs()`
1299: @*/
1300: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1301: {
1302: PetscFunctionBegin;
1304: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1305: PetscFunctionReturn(PETSC_SUCCESS);
1306: }
1308: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
1309: {
1310: PetscFunctionBegin;
1312: VecCheckAssembled(y);
1314: PetscCall(VecSetErrorIfLocked(y, 1));
1315: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1316: if (nv) {
1317: PetscInt zeros = 0;
1319: PetscAssertPointer(alpha, 3);
1320: PetscAssertPointer(x, 4);
1321: for (PetscInt i = 0; i < nv; ++i) {
1325: PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1326: VecCheckSameSize(y, 1, x[i], 4);
1327: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1328: VecCheckAssembled(x[i]);
1329: PetscCall(VecLockReadPush(x[i]));
1330: zeros += alpha[i] == (PetscScalar)0.0;
1331: }
1333: if (zeros < nv) {
1334: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1335: VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1336: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1337: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1338: }
1340: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1341: }
1342: PetscFunctionReturn(PETSC_SUCCESS);
1343: }
1345: /*@
1346: VecMAXPY - Computes `y = y + sum alpha[i] x[i]`
1348: Logically Collective
1350: Input Parameters:
1351: + nv - number of scalars and `x` vectors
1352: . alpha - array of scalars
1353: . y - one vector
1354: - x - array of vectors
1356: Level: intermediate
1358: Notes:
1359: `y` cannot be any of the `x` vectors
1361: The implementation may use BLAS 2 operations when the vectors `y` have been obtained with `VecDuplicateVecs()`
1363: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`, `VecDuplicateVecs()`
1364: @*/
1365: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1366: {
1367: PetscFunctionBegin;
1368: PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1369: PetscFunctionReturn(PETSC_SUCCESS);
1370: }
1372: /*@
1373: VecMAXPBY - Computes `y = beta y + sum alpha[i] x[i]`
1375: Logically Collective
1377: Input Parameters:
1378: + nv - number of scalars and `x` vectors
1379: . alpha - array of scalars
1380: . beta - scalar
1381: . y - one vector
1382: - x - array of vectors
1384: Level: intermediate
1386: Note:
1387: `y` cannot be any of the `x` vectors.
1389: Developer Notes:
1390: This is a convenience routine, but implementations might be able to optimize it, for example, when `beta` is zero.
1392: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1393: @*/
1394: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
1395: {
1396: PetscFunctionBegin;
1398: VecCheckAssembled(y);
1400: PetscCall(VecSetErrorIfLocked(y, 1));
1401: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1404: if (y->ops->maxpby) {
1405: PetscInt zeros = 0;
1407: if (nv) {
1408: PetscAssertPointer(alpha, 3);
1409: PetscAssertPointer(x, 5);
1410: }
1412: for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
1416: PetscCheckSameTypeAndComm(y, 1, x[i], 5);
1417: VecCheckSameSize(y, 1, x[i], 5);
1418: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1419: VecCheckAssembled(x[i]);
1420: PetscCall(VecLockReadPush(x[i]));
1421: zeros += alpha[i] == (PetscScalar)0.0;
1422: }
1424: if (zeros < nv) { // has nonzero alpha
1425: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1426: PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1427: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1428: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1429: } else {
1430: PetscCall(VecScale(y, beta));
1431: }
1433: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1434: } else { // no maxpby
1435: if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1436: else PetscCall(VecScale(y, beta));
1437: PetscCall(VecMAXPY(y, nv, alpha, x));
1438: }
1439: PetscFunctionReturn(PETSC_SUCCESS);
1440: }
1442: /*@
1443: VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
1444: in the order they appear in the array. The concatenated vector resides on the same
1445: communicator and is the same type as the source vectors.
1447: Collective
1449: Input Parameters:
1450: + nx - number of vectors to be concatenated
1451: - X - array containing the vectors to be concatenated in the order of concatenation
1453: Output Parameters:
1454: + Y - concatenated vector
1455: - x_is - array of index sets corresponding to the concatenated components of `Y` (pass `NULL` if not needed)
1457: Level: advanced
1459: Notes:
1460: Concatenation is similar to the functionality of a `VECNEST` object; they both represent combination of
1461: different vector spaces. However, concatenated vectors do not store any information about their
1462: sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1463: manipulation of data in the concatenated vector that corresponds to the original components at creation.
1465: This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1466: has to operate on combined vector spaces and cannot utilize `VECNEST` objects due to incompatibility with
1467: bound projections.
1469: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1470: @*/
1471: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1472: {
1473: MPI_Comm comm;
1474: VecType vec_type;
1475: Vec Ytmp, Xtmp;
1476: IS *is_tmp;
1477: PetscInt i, shift = 0, Xnl, Xng, Xbegin;
1479: PetscFunctionBegin;
1483: PetscAssertPointer(Y, 3);
1485: if ((*X)->ops->concatenate) {
1486: /* use the dedicated concatenation function if available */
1487: PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1488: } else {
1489: /* loop over vectors and start creating IS */
1490: comm = PetscObjectComm((PetscObject)*X);
1491: PetscCall(VecGetType(*X, &vec_type));
1492: PetscCall(PetscMalloc1(nx, &is_tmp));
1493: for (i = 0; i < nx; i++) {
1494: PetscCall(VecGetSize(X[i], &Xng));
1495: PetscCall(VecGetLocalSize(X[i], &Xnl));
1496: PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1497: PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1498: shift += Xng;
1499: }
1500: /* create the concatenated vector */
1501: PetscCall(VecCreate(comm, &Ytmp));
1502: PetscCall(VecSetType(Ytmp, vec_type));
1503: PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1504: PetscCall(VecSetUp(Ytmp));
1505: /* copy data from X array to Y and return */
1506: for (i = 0; i < nx; i++) {
1507: PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1508: PetscCall(VecCopy(X[i], Xtmp));
1509: PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1510: }
1511: *Y = Ytmp;
1512: if (x_is) {
1513: *x_is = is_tmp;
1514: } else {
1515: for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1516: PetscCall(PetscFree(is_tmp));
1517: }
1518: }
1519: PetscFunctionReturn(PETSC_SUCCESS);
1520: }
1522: /* A helper function for VecGetSubVector to check if we can implement it with no-copy (i.e. the subvector shares
1523: memory with the original vector), and the block size of the subvector.
1525: Input Parameters:
1526: + X - the original vector
1527: - is - the index set of the subvector
1529: Output Parameters:
1530: + contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
1531: . start - start of contiguous block, as an offset from the start of the ownership range of the original vector
1532: - blocksize - the block size of the subvector
1534: */
1535: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1536: {
1537: PetscInt gstart, gend, lstart;
1538: PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1539: PetscInt n, N, ibs, vbs, bs = 1;
1541: PetscFunctionBegin;
1542: PetscCall(ISGetLocalSize(is, &n));
1543: PetscCall(ISGetSize(is, &N));
1544: PetscCall(ISGetBlockSize(is, &ibs));
1545: PetscCall(VecGetBlockSize(X, &vbs));
1546: PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1547: PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1548: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1549: if (ibs > 1) {
1550: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1551: bs = ibs;
1552: } else {
1553: if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1554: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1555: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1556: }
1558: *contig = red[0];
1559: *start = lstart;
1560: *blocksize = bs;
1561: PetscFunctionReturn(PETSC_SUCCESS);
1562: }
1564: /* A helper function for VecGetSubVector, to be used when we have to build a standalone subvector through VecScatter
1566: Input Parameters:
1567: + X - the original vector
1568: . is - the index set of the subvector
1569: - bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()
1571: Output Parameter:
1572: . Z - the subvector, which will compose the VecScatter context on output
1573: */
1574: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1575: {
1576: PetscInt n, N;
1577: VecScatter vscat;
1578: Vec Y;
1580: PetscFunctionBegin;
1581: PetscCall(ISGetLocalSize(is, &n));
1582: PetscCall(ISGetSize(is, &N));
1583: PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1584: PetscCall(VecSetSizes(Y, n, N));
1585: PetscCall(VecSetBlockSize(Y, bs));
1586: PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1587: PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1588: PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1589: PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1590: PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1591: PetscCall(VecScatterDestroy(&vscat));
1592: *Z = Y;
1593: PetscFunctionReturn(PETSC_SUCCESS);
1594: }
1596: /*@
1597: VecGetSubVector - Gets a vector representing part of another vector
1599: Collective
1601: Input Parameters:
1602: + X - vector from which to extract a subvector
1603: - is - index set representing portion of `X` to extract
1605: Output Parameter:
1606: . Y - subvector corresponding to `is`
1608: Level: advanced
1610: Notes:
1611: The subvector `Y` should be returned with `VecRestoreSubVector()`.
1612: `X` and `is` must be defined on the same communicator
1614: Changes to the subvector will be reflected in the `X` vector on the call to `VecRestoreSubVector()`.
1616: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1617: modifying the subvector. Other non-overlapping subvectors can still be obtained from `X` using this function.
1619: 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`.
1621: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1622: @*/
1623: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1624: {
1625: Vec Z;
1627: PetscFunctionBegin;
1630: PetscCheckSameComm(X, 1, is, 2);
1631: PetscAssertPointer(Y, 3);
1632: if (X->ops->getsubvector) {
1633: PetscUseTypeMethod(X, getsubvector, is, &Z);
1634: } else { /* Default implementation currently does no caching */
1635: PetscBool contig;
1636: PetscInt n, N, start, bs;
1638: PetscCall(ISGetLocalSize(is, &n));
1639: PetscCall(ISGetSize(is, &N));
1640: PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1641: if (contig) { /* We can do a no-copy implementation */
1642: const PetscScalar *x;
1643: PetscInt state = 0;
1644: PetscBool isstd, iscuda, iship;
1646: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1647: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1648: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1649: if (iscuda) {
1650: #if PetscDefined(HAVE_CUDA)
1651: const PetscScalar *x_d;
1652: PetscMPIInt size;
1653: PetscOffloadMask flg;
1655: PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1656: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1657: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1658: if (x) x += start;
1659: if (x_d) x_d += start;
1660: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1661: if (size == 1) {
1662: PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1663: } else {
1664: PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1665: }
1666: Z->offloadmask = flg;
1667: #endif
1668: } else if (iship) {
1669: #if PetscDefined(HAVE_HIP)
1670: const PetscScalar *x_d;
1671: PetscMPIInt size;
1672: PetscOffloadMask flg;
1674: PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1675: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1676: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1677: if (x) x += start;
1678: if (x_d) x_d += start;
1679: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1680: if (size == 1) {
1681: PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1682: } else {
1683: PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1684: }
1685: Z->offloadmask = flg;
1686: #endif
1687: } else if (isstd) {
1688: PetscMPIInt size;
1690: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1691: PetscCall(VecGetArrayRead(X, &x));
1692: if (x) x += start;
1693: if (size == 1) {
1694: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1695: } else {
1696: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1697: }
1698: PetscCall(VecRestoreArrayRead(X, &x));
1699: } else { /* default implementation: use place array */
1700: PetscCall(VecGetArrayRead(X, &x));
1701: PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1702: PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1703: PetscCall(VecSetSizes(Z, n, N));
1704: PetscCall(VecSetBlockSize(Z, bs));
1705: PetscCall(VecPlaceArray(Z, PetscSafePointerPlusOffset(x, start)));
1706: PetscCall(VecRestoreArrayRead(X, &x));
1707: }
1709: /* this is relevant only in debug mode */
1710: PetscCall(VecLockGet(X, &state));
1711: if (state) PetscCall(VecLockReadPush(Z));
1712: Z->ops->placearray = NULL;
1713: Z->ops->replacearray = NULL;
1714: } else { /* Have to create a scatter and do a copy */
1715: PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1716: }
1717: }
1718: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1719: if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1720: PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1721: *Y = Z;
1722: PetscFunctionReturn(PETSC_SUCCESS);
1723: }
1725: /*@
1726: VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`
1728: Collective
1730: Input Parameters:
1731: + X - vector from which subvector was obtained
1732: . is - index set representing the subset of `X`
1733: - Y - subvector being restored
1735: Level: advanced
1737: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1738: @*/
1739: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1740: {
1741: PETSC_UNUSED PetscObjectState dummystate = 0;
1742: PetscBool unchanged;
1744: PetscFunctionBegin;
1747: PetscCheckSameComm(X, 1, is, 2);
1748: PetscAssertPointer(Y, 3);
1751: if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1752: else {
1753: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1754: if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1755: VecScatter scatter;
1756: PetscInt state;
1758: PetscCall(VecLockGet(X, &state));
1759: PetscCheck(state == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec X is locked for read-only or read/write access");
1761: PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1762: if (scatter) {
1763: PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1764: PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1765: } else {
1766: PetscBool iscuda, iship;
1767: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1768: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1770: if (iscuda) {
1771: #if PetscDefined(HAVE_CUDA)
1772: PetscOffloadMask ymask = (*Y)->offloadmask;
1774: /* The offloadmask of X dictates where to move memory
1775: If X GPU data is valid, then move Y data on GPU if needed
1776: Otherwise, move back to the CPU */
1777: switch (X->offloadmask) {
1778: case PETSC_OFFLOAD_BOTH:
1779: if (ymask == PETSC_OFFLOAD_CPU) {
1780: PetscCall(VecCUDAResetArray(*Y));
1781: } else if (ymask == PETSC_OFFLOAD_GPU) {
1782: X->offloadmask = PETSC_OFFLOAD_GPU;
1783: }
1784: break;
1785: case PETSC_OFFLOAD_GPU:
1786: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1787: break;
1788: case PETSC_OFFLOAD_CPU:
1789: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1790: break;
1791: case PETSC_OFFLOAD_UNALLOCATED:
1792: case PETSC_OFFLOAD_KOKKOS:
1793: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1794: }
1795: #endif
1796: } else if (iship) {
1797: #if PetscDefined(HAVE_HIP)
1798: PetscOffloadMask ymask = (*Y)->offloadmask;
1800: /* The offloadmask of X dictates where to move memory
1801: If X GPU data is valid, then move Y data on GPU if needed
1802: Otherwise, move back to the CPU */
1803: switch (X->offloadmask) {
1804: case PETSC_OFFLOAD_BOTH:
1805: if (ymask == PETSC_OFFLOAD_CPU) {
1806: PetscCall(VecHIPResetArray(*Y));
1807: } else if (ymask == PETSC_OFFLOAD_GPU) {
1808: X->offloadmask = PETSC_OFFLOAD_GPU;
1809: }
1810: break;
1811: case PETSC_OFFLOAD_GPU:
1812: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1813: break;
1814: case PETSC_OFFLOAD_CPU:
1815: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1816: break;
1817: case PETSC_OFFLOAD_UNALLOCATED:
1818: case PETSC_OFFLOAD_KOKKOS:
1819: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1820: }
1821: #endif
1822: } else {
1823: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1824: PetscCall(VecResetArray(*Y));
1825: }
1826: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1827: }
1828: }
1829: }
1830: PetscCall(VecDestroy(Y));
1831: PetscFunctionReturn(PETSC_SUCCESS);
1832: }
1834: /*@
1835: VecCreateLocalVector - Creates a vector object suitable for use with `VecGetLocalVector()` and friends. You must call `VecDestroy()` when the
1836: vector is no longer needed.
1838: Not Collective.
1840: Input Parameter:
1841: . v - The vector for which the local vector is desired.
1843: Output Parameter:
1844: . w - Upon exit this contains the local vector.
1846: Level: beginner
1848: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1849: @*/
1850: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1851: {
1852: VecType roottype;
1853: PetscInt n;
1855: PetscFunctionBegin;
1857: PetscAssertPointer(w, 2);
1858: if (v->ops->createlocalvector) {
1859: PetscUseTypeMethod(v, createlocalvector, w);
1860: PetscFunctionReturn(PETSC_SUCCESS);
1861: }
1862: PetscCall(VecGetRootType_Private(v, &roottype));
1863: PetscCall(VecCreate(PETSC_COMM_SELF, w));
1864: PetscCall(VecGetLocalSize(v, &n));
1865: PetscCall(VecSetSizes(*w, n, n));
1866: PetscCall(VecGetBlockSize(v, &n));
1867: PetscCall(VecSetBlockSize(*w, n));
1868: PetscCall(VecSetType(*w, roottype));
1869: PetscFunctionReturn(PETSC_SUCCESS);
1870: }
1872: /*@
1873: VecGetLocalVectorRead - Maps the local portion of a vector into a
1874: vector.
1876: Not Collective.
1878: Input Parameter:
1879: . v - The vector for which the local vector is desired.
1881: Output Parameter:
1882: . w - Upon exit this contains the local vector.
1884: Level: beginner
1886: Notes:
1887: You must call `VecRestoreLocalVectorRead()` when the local
1888: vector is no longer needed.
1890: This function is similar to `VecGetArrayRead()` which maps the local
1891: portion into a raw pointer. `VecGetLocalVectorRead()` is usually
1892: almost as efficient as `VecGetArrayRead()` but in certain circumstances
1893: `VecGetLocalVectorRead()` can be much more efficient than
1894: `VecGetArrayRead()`. This is because the construction of a contiguous
1895: array representing the vector data required by `VecGetArrayRead()` can
1896: be an expensive operation for certain vector types. For example, for
1897: GPU vectors `VecGetArrayRead()` requires that the data between device
1898: and host is synchronized.
1900: Unlike `VecGetLocalVector()`, this routine is not collective and
1901: preserves cached information.
1903: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1904: @*/
1905: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1906: {
1907: PetscFunctionBegin;
1910: VecCheckSameLocalSize(v, 1, w, 2);
1911: if (v->ops->getlocalvectorread) {
1912: PetscUseTypeMethod(v, getlocalvectorread, w);
1913: } else {
1914: PetscScalar *a;
1916: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1917: PetscCall(VecPlaceArray(w, a));
1918: }
1919: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1920: PetscCall(VecLockReadPush(v));
1921: PetscCall(VecLockReadPush(w));
1922: PetscFunctionReturn(PETSC_SUCCESS);
1923: }
1925: /*@
1926: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1927: previously mapped into a vector using `VecGetLocalVectorRead()`.
1929: Not Collective.
1931: Input Parameters:
1932: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
1933: - w - The vector into which the local portion of `v` was mapped.
1935: Level: beginner
1937: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1938: @*/
1939: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1940: {
1941: PetscFunctionBegin;
1944: if (v->ops->restorelocalvectorread) {
1945: PetscUseTypeMethod(v, restorelocalvectorread, w);
1946: } else {
1947: const PetscScalar *a;
1949: PetscCall(VecGetArrayRead(w, &a));
1950: PetscCall(VecRestoreArrayRead(v, &a));
1951: PetscCall(VecResetArray(w));
1952: }
1953: PetscCall(VecLockReadPop(v));
1954: PetscCall(VecLockReadPop(w));
1955: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1956: PetscFunctionReturn(PETSC_SUCCESS);
1957: }
1959: /*@
1960: VecGetLocalVector - Maps the local portion of a vector into a
1961: vector.
1963: Collective
1965: Input Parameter:
1966: . v - The vector for which the local vector is desired.
1968: Output Parameter:
1969: . w - Upon exit this contains the local vector.
1971: Level: beginner
1973: Notes:
1974: You must call `VecRestoreLocalVector()` when the local
1975: vector is no longer needed.
1977: This function is similar to `VecGetArray()` which maps the local
1978: portion into a raw pointer. `VecGetLocalVector()` is usually about as
1979: efficient as `VecGetArray()` but in certain circumstances
1980: `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1981: This is because the construction of a contiguous array representing
1982: the vector data required by `VecGetArray()` can be an expensive
1983: operation for certain vector types. For example, for GPU vectors
1984: `VecGetArray()` requires that the data between device and host is
1985: synchronized.
1987: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1988: @*/
1989: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1990: {
1991: PetscFunctionBegin;
1994: VecCheckSameLocalSize(v, 1, w, 2);
1995: if (v->ops->getlocalvector) {
1996: PetscUseTypeMethod(v, getlocalvector, w);
1997: } else {
1998: PetscScalar *a;
2000: PetscCall(VecGetArray(v, &a));
2001: PetscCall(VecPlaceArray(w, a));
2002: }
2003: PetscCall(PetscObjectStateIncrease((PetscObject)w));
2004: PetscFunctionReturn(PETSC_SUCCESS);
2005: }
2007: /*@
2008: VecRestoreLocalVector - Unmaps the local portion of a vector
2009: previously mapped into a vector using `VecGetLocalVector()`.
2011: Logically Collective.
2013: Input Parameters:
2014: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
2015: - w - The vector into which the local portion of `v` was mapped.
2017: Level: beginner
2019: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
2020: @*/
2021: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
2022: {
2023: PetscFunctionBegin;
2026: if (v->ops->restorelocalvector) {
2027: PetscUseTypeMethod(v, restorelocalvector, w);
2028: } else {
2029: PetscScalar *a;
2030: PetscCall(VecGetArray(w, &a));
2031: PetscCall(VecRestoreArray(v, &a));
2032: PetscCall(VecResetArray(w));
2033: }
2034: PetscCall(PetscObjectStateIncrease((PetscObject)w));
2035: PetscCall(PetscObjectStateIncrease((PetscObject)v));
2036: PetscFunctionReturn(PETSC_SUCCESS);
2037: }
2039: /*@
2040: VecGetArray - Returns a pointer to a contiguous array that contains this
2041: MPI processes's portion of the vector data
2043: Logically Collective
2045: Input Parameter:
2046: . x - the vector
2048: Output Parameter:
2049: . a - location to put pointer to the array
2051: Level: beginner
2053: Notes:
2054: For the standard PETSc vectors, `VecGetArray()` returns a pointer to the local data array and
2055: does not use any copies. If the underlying vector data is not stored in a contiguous array
2056: this routine will copy the data to a contiguous array and return a pointer to that. You MUST
2057: call `VecRestoreArray()` when you no longer need access to the array.
2059: For vectors that may also have the array data in GPU memory, for example, `VECCUDA`, this call ensures the CPU array has the
2060: most recent array values by copying the data from the GPU memory if needed.
2062: Fortran Note:
2063: .vb
2064: PetscScalar, pointer :: a(:)
2065: .ve
2067: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecPlaceArray()`, `VecGetArray2d()`,
2068: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`, `VecGetArrayAndMemType()`
2069: @*/
2070: PetscErrorCode VecGetArray(Vec x, PetscScalar *a[])
2071: {
2072: PetscFunctionBegin;
2074: PetscCall(VecSetErrorIfLocked(x, 1));
2075: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
2076: PetscUseTypeMethod(x, getarray, a);
2077: } else if (x->petscnative) { /* VECSTANDARD */
2078: *a = *((PetscScalar **)x->data);
2079: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
2080: PetscFunctionReturn(PETSC_SUCCESS);
2081: }
2083: /*@
2084: VecRestoreArray - Restores a vector after `VecGetArray()` has been called and the array is no longer needed
2086: Logically Collective
2088: Input Parameters:
2089: + x - the vector
2090: - a - location of pointer to array obtained from `VecGetArray()`
2092: Level: beginner
2094: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2095: `VecGetArrayPair()`, `VecRestoreArrayPair()`
2096: @*/
2097: PetscErrorCode VecRestoreArray(Vec x, PetscScalar *a[])
2098: {
2099: PetscFunctionBegin;
2101: if (a) PetscAssertPointer(a, 2);
2102: if (x->ops->restorearray) {
2103: PetscUseTypeMethod(x, restorearray, a);
2104: } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
2105: if (a) *a = NULL;
2106: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2107: PetscFunctionReturn(PETSC_SUCCESS);
2108: }
2109: /*@
2110: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
2112: Not Collective
2114: Input Parameter:
2115: . x - the vector
2117: Output Parameter:
2118: . a - the array
2120: Level: beginner
2122: Notes:
2123: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2125: Unlike `VecGetArray()`, preserves cached information like vector norms.
2127: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
2128: implementations may require a copy, but such implementations should cache the contiguous representation so that
2129: only one copy is performed when this routine is called multiple times in sequence.
2131: For vectors that may also have the array data in GPU memory, for example, `VECCUDA`, this call ensures the CPU array has the
2132: most recent array values by copying the data from the GPU memory if needed.
2134: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2135: `VecGetArrayAndMemType()`
2136: @*/
2137: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar *a[])
2138: {
2139: PetscFunctionBegin;
2141: PetscAssertPointer(a, 2);
2142: if (x->ops->getarrayread) {
2143: PetscUseTypeMethod(x, getarrayread, a);
2144: } else if (x->ops->getarray) {
2145: PetscObjectState state;
2147: /* VECNEST, VECCUDA, VECKOKKOS etc */
2148: // x->ops->getarray may bump the object state, but since we know this is a read-only get
2149: // we can just undo that
2150: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2151: PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
2152: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2153: } else if (x->petscnative) {
2154: /* VECSTANDARD */
2155: *a = *((PetscScalar **)x->data);
2156: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
2157: PetscFunctionReturn(PETSC_SUCCESS);
2158: }
2160: /*@
2161: VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`
2163: Not Collective
2165: Input Parameters:
2166: + x - the vector
2167: - a - the array
2169: Level: beginner
2171: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2172: @*/
2173: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar *a[])
2174: {
2175: PetscFunctionBegin;
2177: if (a) PetscAssertPointer(a, 2);
2178: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
2179: /* nothing */
2180: } else if (x->ops->restorearrayread) { /* VECNEST */
2181: PetscUseTypeMethod(x, restorearrayread, a);
2182: } else { /* No one? */
2183: PetscObjectState state;
2185: // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2186: // we can just undo that
2187: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2188: PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2189: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2190: }
2191: if (a) *a = NULL;
2192: PetscFunctionReturn(PETSC_SUCCESS);
2193: }
2195: /*@
2196: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
2197: MPI processes's portion of the vector data.
2199: Logically Collective
2201: Input Parameter:
2202: . x - the vector
2204: Output Parameter:
2205: . a - location to put pointer to the array
2207: Level: intermediate
2209: Note:
2210: The values in this array are NOT valid, the caller of this routine is responsible for putting
2211: values into the array; any values it does not set will be invalid.
2213: The array must be returned using a matching call to `VecRestoreArrayWrite()`.
2215: For vectors associated with GPUs, the host and device vectors are not synchronized before
2216: giving access. If you need correct values in the array use `VecGetArray()`
2218: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecPlaceArray()`, `VecGetArray2d()`,
2219: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`, `VecGetArrayAndMemType()`
2220: @*/
2221: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar *a[])
2222: {
2223: PetscFunctionBegin;
2225: PetscAssertPointer(a, 2);
2226: PetscCall(VecSetErrorIfLocked(x, 1));
2227: if (x->ops->getarraywrite) {
2228: PetscUseTypeMethod(x, getarraywrite, a);
2229: } else {
2230: PetscCall(VecGetArray(x, a));
2231: }
2232: PetscFunctionReturn(PETSC_SUCCESS);
2233: }
2235: /*@
2236: VecRestoreArrayWrite - Restores a vector after `VecGetArrayWrite()` has been called.
2238: Logically Collective
2240: Input Parameters:
2241: + x - the vector
2242: - a - location of pointer to array obtained from `VecGetArray()`
2244: Level: beginner
2246: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2247: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2248: @*/
2249: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar *a[])
2250: {
2251: PetscFunctionBegin;
2253: if (a) PetscAssertPointer(a, 2);
2254: if (x->ops->restorearraywrite) {
2255: PetscUseTypeMethod(x, restorearraywrite, a);
2256: } else if (x->ops->restorearray) {
2257: PetscUseTypeMethod(x, restorearray, a);
2258: }
2259: if (a) *a = NULL;
2260: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2261: PetscFunctionReturn(PETSC_SUCCESS);
2262: }
2264: /*@
2265: VecGetArrays - Returns a pointer to the arrays in a set of vectors
2266: that were created by a call to `VecDuplicateVecs()`.
2268: Logically Collective; No Fortran Support
2270: Input Parameters:
2271: + x - the vectors
2272: - n - the number of vectors
2274: Output Parameter:
2275: . a - location to put pointer to the array
2277: Level: intermediate
2279: Note:
2280: You MUST call `VecRestoreArrays()` when you no longer need access to the arrays.
2282: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2283: @*/
2284: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2285: {
2286: PetscInt i;
2287: PetscScalar **q;
2289: PetscFunctionBegin;
2290: PetscAssertPointer(x, 1);
2292: PetscAssertPointer(a, 3);
2293: PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2294: PetscCall(PetscMalloc1(n, &q));
2295: for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2296: *a = q;
2297: PetscFunctionReturn(PETSC_SUCCESS);
2298: }
2300: /*@
2301: VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2302: has been called.
2304: Logically Collective; No Fortran Support
2306: Input Parameters:
2307: + x - the vector
2308: . n - the number of vectors
2309: - a - location of pointer to arrays obtained from `VecGetArrays()`
2311: Notes:
2312: For regular PETSc vectors this routine does not involve any copies. For
2313: any special vectors that do not store local vector data in a contiguous
2314: array, this routine will copy the data back into the underlying
2315: vector data structure from the arrays obtained with `VecGetArrays()`.
2317: Level: intermediate
2319: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2320: @*/
2321: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2322: {
2323: PetscScalar **q = *a;
2325: PetscFunctionBegin;
2326: PetscAssertPointer(x, 1);
2328: PetscAssertPointer(a, 3);
2330: for (PetscInt i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2331: PetscCall(PetscFree(q));
2332: PetscFunctionReturn(PETSC_SUCCESS);
2333: }
2335: /*@
2336: VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g.,
2337: `VECCUDA`), the returned pointer will be a device pointer to the device memory that contains
2338: this MPI processes's portion of the vector data.
2340: Logically Collective; No Fortran Support
2342: Input Parameter:
2343: . x - the vector
2345: Output Parameters:
2346: + a - location to put pointer to the array
2347: - mtype - memory type of the array
2349: Level: beginner
2351: Note:
2352: Device data is guaranteed to have the latest value. Otherwise, when this is a host vector
2353: (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host
2354: pointer.
2356: For `VECKOKKOS`, if Kokkos is configured without device (e.g., use serial or openmp), per
2357: this function, the vector works like `VECSEQ`/`VECMPI`; otherwise, it works like `VECCUDA` or
2358: `VECHIP` etc.
2360: Use `VecRestoreArrayAndMemType()` when the array access is no longer needed.
2362: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`,
2363: `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2364: @*/
2365: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar *a[], PetscMemType *mtype)
2366: {
2367: PetscFunctionBegin;
2370: if (a) PetscAssertPointer(a, 2);
2371: if (mtype) PetscAssertPointer(mtype, 3);
2372: PetscCall(VecSetErrorIfLocked(x, 1));
2373: if (x->ops->getarrayandmemtype) {
2374: /* VECCUDA, VECKOKKOS etc */
2375: PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2376: } else {
2377: /* VECSTANDARD, VECNEST, VECVIENNACL */
2378: PetscCall(VecGetArray(x, a));
2379: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2380: }
2381: PetscFunctionReturn(PETSC_SUCCESS);
2382: }
2384: /*@
2385: VecRestoreArrayAndMemType - Restores a vector after `VecGetArrayAndMemType()` has been called.
2387: Logically Collective; No Fortran Support
2389: Input Parameters:
2390: + x - the vector
2391: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`
2393: Level: beginner
2395: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`,
2396: `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2397: @*/
2398: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar *a[])
2399: {
2400: PetscFunctionBegin;
2403: if (a) PetscAssertPointer(a, 2);
2404: if (x->ops->restorearrayandmemtype) {
2405: /* VECCUDA, VECKOKKOS etc */
2406: PetscUseTypeMethod(x, restorearrayandmemtype, a);
2407: } else {
2408: /* VECNEST, VECVIENNACL */
2409: PetscCall(VecRestoreArray(x, a));
2410: } /* VECSTANDARD does nothing */
2411: if (a) *a = NULL;
2412: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2413: PetscFunctionReturn(PETSC_SUCCESS);
2414: }
2416: /*@
2417: VecGetArrayReadAndMemType - Like `VecGetArrayRead()`, but if the input vector is a device vector, it will return a read-only device pointer.
2418: The returned pointer is guaranteed to point to up-to-date data. For host vectors, it functions as `VecGetArrayRead()`.
2420: Not Collective; No Fortran Support
2422: Input Parameter:
2423: . x - the vector
2425: Output Parameters:
2426: + a - the array
2427: - mtype - memory type of the array
2429: Level: beginner
2431: Notes:
2432: The array must be returned using a matching call to `VecRestoreArrayReadAndMemType()`.
2434: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2435: @*/
2436: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar *a[], PetscMemType *mtype)
2437: {
2438: PetscFunctionBegin;
2441: PetscAssertPointer(a, 2);
2442: if (mtype) PetscAssertPointer(mtype, 3);
2443: if (x->ops->getarrayreadandmemtype) {
2444: /* VECCUDA/VECHIP though they are also petscnative */
2445: PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2446: } else if (x->ops->getarrayandmemtype) {
2447: /* VECKOKKOS */
2448: PetscObjectState state;
2450: // see VecGetArrayRead() for why
2451: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2452: PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2453: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2454: } else {
2455: PetscCall(VecGetArrayRead(x, a));
2456: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2457: }
2458: PetscFunctionReturn(PETSC_SUCCESS);
2459: }
2461: /*@
2462: VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`
2464: Not Collective; No Fortran Support
2466: Input Parameters:
2467: + x - the vector
2468: - a - the array
2470: Level: beginner
2472: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2473: @*/
2474: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar *a[])
2475: {
2476: PetscFunctionBegin;
2479: if (a) PetscAssertPointer(a, 2);
2480: if (x->ops->restorearrayreadandmemtype) {
2481: /* VECCUDA/VECHIP */
2482: PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2483: } else if (!x->petscnative) {
2484: /* VECNEST */
2485: PetscCall(VecRestoreArrayRead(x, a));
2486: }
2487: if (a) *a = NULL;
2488: PetscFunctionReturn(PETSC_SUCCESS);
2489: }
2491: /*@
2492: VecGetArrayWriteAndMemType - Like `VecGetArrayWrite()`, but if this is a device vector it will always return
2493: a device pointer to the device memory that contains this processor's portion of the vector data.
2495: Logically Collective; No Fortran Support
2497: Input Parameter:
2498: . x - the vector
2500: Output Parameters:
2501: + a - the array
2502: - mtype - memory type of the array
2504: Level: beginner
2506: Note:
2507: The array must be returned using a matching call to `VecRestoreArrayWriteAndMemType()`, where it will label the device memory as most recent.
2509: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2510: @*/
2511: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar *a[], PetscMemType *mtype)
2512: {
2513: PetscFunctionBegin;
2516: PetscCall(VecSetErrorIfLocked(x, 1));
2517: PetscAssertPointer(a, 2);
2518: if (mtype) PetscAssertPointer(mtype, 3);
2519: if (x->ops->getarraywriteandmemtype) {
2520: /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2521: PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2522: } else if (x->ops->getarrayandmemtype) {
2523: PetscCall(VecGetArrayAndMemType(x, a, mtype));
2524: } else {
2525: /* VECNEST, VECVIENNACL */
2526: PetscCall(VecGetArrayWrite(x, a));
2527: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2528: }
2529: PetscFunctionReturn(PETSC_SUCCESS);
2530: }
2532: /*@
2533: VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`
2535: Logically Collective; No Fortran Support
2537: Input Parameters:
2538: + x - the vector
2539: - a - the array
2541: Level: beginner
2543: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2544: @*/
2545: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar *a[])
2546: {
2547: PetscFunctionBegin;
2550: PetscCall(VecSetErrorIfLocked(x, 1));
2551: if (a) PetscAssertPointer(a, 2);
2552: if (x->ops->restorearraywriteandmemtype) {
2553: /* VECCUDA/VECHIP */
2554: PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2555: PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2556: } else if (x->ops->restorearrayandmemtype) {
2557: PetscCall(VecRestoreArrayAndMemType(x, a));
2558: } else {
2559: PetscCall(VecRestoreArray(x, a));
2560: }
2561: if (a) *a = NULL;
2562: PetscFunctionReturn(PETSC_SUCCESS);
2563: }
2565: /*@
2566: VecPlaceArray - Allows one to replace the array in a vector with an
2567: array provided by the user. This is useful to avoid copying an array
2568: into a vector.
2570: Logically Collective
2572: Input Parameters:
2573: + vec - the vector
2574: - array - the array
2576: Level: developer
2578: Notes:
2579: Adding `const` to `array` was an oversight, as subsequent operations on `vec` would
2580: likely modify the data in `array`. However, we have kept it to avoid breaking APIs.
2582: Use `VecReplaceArray()` instead to permanently replace the array
2584: You can return to the original array with a call to `VecResetArray()`. `vec` does not take
2585: ownership of `array` in any way.
2587: The user must free `array` themselves but be careful not to
2588: do so before the vector has either been destroyed, had its original array restored with
2589: `VecResetArray()` or permanently replaced with `VecReplaceArray()`.
2591: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2592: @*/
2593: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2594: {
2595: PetscFunctionBegin;
2598: if (array) PetscAssertPointer(array, 2);
2599: PetscUseTypeMethod(vec, placearray, array);
2600: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2601: PetscFunctionReturn(PETSC_SUCCESS);
2602: }
2604: /*@
2605: VecReplaceArray - Allows one to replace the array in a vector with an
2606: array provided by the user. This is useful to avoid copying an array
2607: into a vector.
2609: Logically Collective; No Fortran Support
2611: Input Parameters:
2612: + vec - the vector
2613: - array - the array
2615: Level: developer
2617: Notes:
2618: Adding `const` to `array` was an oversight, as subsequent operations on `vec` would
2619: likely modify the data in `array`. However, we have kept it to avoid breaking APIs.
2621: This permanently replaces the array and frees the memory associated
2622: with the old array. Use `VecPlaceArray()` to temporarily replace the array.
2624: The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
2625: freed by the user. It will be freed when the vector is destroyed.
2627: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2628: @*/
2629: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2630: {
2631: PetscFunctionBegin;
2634: PetscUseTypeMethod(vec, replacearray, array);
2635: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2636: PetscFunctionReturn(PETSC_SUCCESS);
2637: }
2639: /*@
2640: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2641: processor's portion of the vector data. You MUST call `VecRestoreArray2d()`
2642: when you no longer need access to the array.
2644: Logically Collective
2646: Input Parameters:
2647: + x - the vector
2648: . m - first dimension of two dimensional array
2649: . n - second dimension of two dimensional array
2650: . mstart - first index you will use in first coordinate direction (often 0)
2651: - nstart - first index in the second coordinate direction (often 0)
2653: Output Parameter:
2654: . a - location to put pointer to the array
2656: Level: developer
2658: Notes:
2659: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2660: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2661: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2662: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2664: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2666: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2667: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2668: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2669: @*/
2670: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2671: {
2672: PetscInt i, N;
2673: PetscScalar *aa;
2675: PetscFunctionBegin;
2677: PetscAssertPointer(a, 6);
2679: PetscCall(VecGetLocalSize(x, &N));
2680: 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);
2681: PetscCall(VecGetArray(x, &aa));
2683: PetscCall(PetscMalloc1(m, a));
2684: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2685: *a -= mstart;
2686: PetscFunctionReturn(PETSC_SUCCESS);
2687: }
2689: /*@
2690: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2691: processor's portion of the vector data. You MUST call `VecRestoreArray2dWrite()`
2692: when you no longer need access to the array.
2694: Logically Collective
2696: Input Parameters:
2697: + x - the vector
2698: . m - first dimension of two dimensional array
2699: . n - second dimension of two dimensional array
2700: . mstart - first index you will use in first coordinate direction (often 0)
2701: - nstart - first index in the second coordinate direction (often 0)
2703: Output Parameter:
2704: . a - location to put pointer to the array
2706: Level: developer
2708: Notes:
2709: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2710: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2711: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2712: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2714: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2716: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2717: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2718: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2719: @*/
2720: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2721: {
2722: PetscInt i, N;
2723: PetscScalar *aa;
2725: PetscFunctionBegin;
2727: PetscAssertPointer(a, 6);
2729: PetscCall(VecGetLocalSize(x, &N));
2730: 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);
2731: PetscCall(VecGetArrayWrite(x, &aa));
2733: PetscCall(PetscMalloc1(m, a));
2734: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2735: *a -= mstart;
2736: PetscFunctionReturn(PETSC_SUCCESS);
2737: }
2739: /*@
2740: VecRestoreArray2d - Restores a vector after `VecGetArray2d()` 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 VecRestoreArray2d(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(VecRestoreArray(x, NULL));
2777: *a = NULL;
2778: PetscFunctionReturn(PETSC_SUCCESS);
2779: }
2781: /*@
2782: VecRestoreArray2dWrite - Restores a vector after `VecGetArray2dWrite()` has been called.
2784: Logically Collective
2786: Input Parameters:
2787: + x - the vector
2788: . m - first dimension of two dimensional array
2789: . n - second dimension of the two dimensional array
2790: . mstart - first index you will use in first coordinate direction (often 0)
2791: . nstart - first index in the second coordinate direction (often 0)
2792: - a - location of pointer to array obtained from `VecGetArray2d()`
2794: Level: developer
2796: Notes:
2797: For regular PETSc vectors this routine does not involve any copies. For
2798: any special vectors that do not store local vector data in a contiguous
2799: array, this routine will copy the data back into the underlying
2800: vector data structure from the array obtained with `VecGetArray()`.
2802: This routine actually zeros out the `a` pointer.
2804: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2805: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
2806: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2807: @*/
2808: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2809: {
2810: void *dummy;
2812: PetscFunctionBegin;
2814: PetscAssertPointer(a, 6);
2816: dummy = (void *)(*a + mstart);
2817: PetscCall(PetscFree(dummy));
2818: PetscCall(VecRestoreArrayWrite(x, NULL));
2819: PetscFunctionReturn(PETSC_SUCCESS);
2820: }
2822: /*@
2823: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2824: processor's portion of the vector data. You MUST call `VecRestoreArray1d()`
2825: when you no longer need access to the array.
2827: Logically Collective
2829: Input Parameters:
2830: + x - the vector
2831: . m - first dimension of two dimensional array
2832: - mstart - first index you will use in first coordinate direction (often 0)
2834: Output Parameter:
2835: . a - location to put pointer to the array
2837: Level: developer
2839: Notes:
2840: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2841: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2842: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2844: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2846: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2847: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2848: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2849: @*/
2850: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2851: {
2852: PetscInt N;
2854: PetscFunctionBegin;
2856: PetscAssertPointer(a, 4);
2858: PetscCall(VecGetLocalSize(x, &N));
2859: 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);
2860: PetscCall(VecGetArray(x, a));
2861: *a -= mstart;
2862: PetscFunctionReturn(PETSC_SUCCESS);
2863: }
2865: /*@
2866: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
2867: processor's portion of the vector data. You MUST call `VecRestoreArray1dWrite()`
2868: when you no longer need access to the array.
2870: Logically Collective
2872: Input Parameters:
2873: + x - the vector
2874: . m - first dimension of two dimensional array
2875: - mstart - first index you will use in first coordinate direction (often 0)
2877: Output Parameter:
2878: . a - location to put pointer to the array
2880: Level: developer
2882: Notes:
2883: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2884: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2885: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2887: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2889: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2890: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2891: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2892: @*/
2893: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2894: {
2895: PetscInt N;
2897: PetscFunctionBegin;
2899: PetscAssertPointer(a, 4);
2901: PetscCall(VecGetLocalSize(x, &N));
2902: 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);
2903: PetscCall(VecGetArrayWrite(x, a));
2904: *a -= mstart;
2905: PetscFunctionReturn(PETSC_SUCCESS);
2906: }
2908: /*@
2909: VecRestoreArray1d - Restores a vector after `VecGetArray1d()` has been called.
2911: Logically Collective
2913: Input Parameters:
2914: + x - the vector
2915: . m - first dimension of two dimensional array
2916: . mstart - first index you will use in first coordinate direction (often 0)
2917: - a - location of pointer to array obtained from `VecGetArray1d()`
2919: Level: developer
2921: Notes:
2922: For regular PETSc vectors this routine does not involve any copies. For
2923: any special vectors that do not store local vector data in a contiguous
2924: array, this routine will copy the data back into the underlying
2925: vector data structure from the array obtained with `VecGetArray1d()`.
2927: This routine actually zeros out the `a` pointer.
2929: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2930: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
2931: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2932: @*/
2933: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2934: {
2935: PetscFunctionBegin;
2938: PetscCall(VecRestoreArray(x, NULL));
2939: *a = NULL;
2940: PetscFunctionReturn(PETSC_SUCCESS);
2941: }
2943: /*@
2944: VecRestoreArray1dWrite - Restores a vector after `VecGetArray1dWrite()` has been called.
2946: Logically Collective
2948: Input Parameters:
2949: + x - the vector
2950: . m - first dimension of two dimensional array
2951: . mstart - first index you will use in first coordinate direction (often 0)
2952: - a - location of pointer to array obtained from `VecGetArray1d()`
2954: Level: developer
2956: Notes:
2957: For regular PETSc vectors this routine does not involve any copies. For
2958: any special vectors that do not store local vector data in a contiguous
2959: array, this routine will copy the data back into the underlying
2960: vector data structure from the array obtained with `VecGetArray1d()`.
2962: This routine actually zeros out the `a` pointer.
2964: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2965: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
2966: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2967: @*/
2968: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2969: {
2970: PetscFunctionBegin;
2973: PetscCall(VecRestoreArrayWrite(x, NULL));
2974: *a = NULL;
2975: PetscFunctionReturn(PETSC_SUCCESS);
2976: }
2978: /*@
2979: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
2980: processor's portion of the vector data. You MUST call `VecRestoreArray3d()`
2981: when you no longer need access to the array.
2983: Logically Collective
2985: Input Parameters:
2986: + x - the vector
2987: . m - first dimension of three dimensional array
2988: . n - second dimension of three dimensional array
2989: . p - third dimension of three dimensional array
2990: . mstart - first index you will use in first coordinate direction (often 0)
2991: . nstart - first index in the second coordinate direction (often 0)
2992: - pstart - first index in the third coordinate direction (often 0)
2994: Output Parameter:
2995: . a - location to put pointer to the array
2997: Level: developer
2999: Notes:
3000: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3001: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3002: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3003: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3005: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3007: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3008: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
3009: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3010: @*/
3011: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3012: {
3013: PetscInt i, N, j;
3014: PetscScalar *aa, **b;
3016: PetscFunctionBegin;
3018: PetscAssertPointer(a, 8);
3020: PetscCall(VecGetLocalSize(x, &N));
3021: 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);
3022: PetscCall(VecGetArray(x, &aa));
3024: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3025: b = (PetscScalar **)((*a) + m);
3026: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3027: for (i = 0; i < m; i++)
3028: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset(aa, i * n * p + j * p - pstart);
3029: *a -= mstart;
3030: PetscFunctionReturn(PETSC_SUCCESS);
3031: }
3033: /*@
3034: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
3035: processor's portion of the vector data. You MUST call `VecRestoreArray3dWrite()`
3036: when you no longer need access to the array.
3038: Logically Collective
3040: Input Parameters:
3041: + x - the vector
3042: . m - first dimension of three dimensional array
3043: . n - second dimension of three dimensional array
3044: . p - third dimension of three dimensional array
3045: . mstart - first index you will use in first coordinate direction (often 0)
3046: . nstart - first index in the second coordinate direction (often 0)
3047: - pstart - first index in the third coordinate direction (often 0)
3049: Output Parameter:
3050: . a - location to put pointer to the array
3052: Level: developer
3054: Notes:
3055: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3056: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3057: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3058: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3060: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3062: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3063: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3064: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3065: @*/
3066: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3067: {
3068: PetscInt i, N, j;
3069: PetscScalar *aa, **b;
3071: PetscFunctionBegin;
3073: PetscAssertPointer(a, 8);
3075: PetscCall(VecGetLocalSize(x, &N));
3076: 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);
3077: PetscCall(VecGetArrayWrite(x, &aa));
3079: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3080: b = (PetscScalar **)((*a) + m);
3081: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3082: for (i = 0; i < m; i++)
3083: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3085: *a -= mstart;
3086: PetscFunctionReturn(PETSC_SUCCESS);
3087: }
3089: /*@
3090: VecRestoreArray3d - Restores a vector after `VecGetArray3d()` has been called.
3092: Logically Collective
3094: Input Parameters:
3095: + x - the vector
3096: . m - first dimension of three dimensional array
3097: . n - second dimension of the three dimensional array
3098: . p - third dimension of the three dimensional array
3099: . mstart - first index you will use in first coordinate direction (often 0)
3100: . nstart - first index in the second coordinate direction (often 0)
3101: . pstart - first index in the third coordinate direction (often 0)
3102: - a - location of pointer to array obtained from VecGetArray3d()
3104: Level: developer
3106: Notes:
3107: For regular PETSc vectors this routine does not involve any copies. For
3108: any special vectors that do not store local vector data in a contiguous
3109: array, this routine will copy the data back into the underlying
3110: vector data structure from the array obtained with `VecGetArray()`.
3112: This routine actually zeros out the `a` pointer.
3114: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3115: `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3116: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3117: @*/
3118: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3119: {
3120: void *dummy;
3122: PetscFunctionBegin;
3124: PetscAssertPointer(a, 8);
3126: dummy = (void *)(*a + mstart);
3127: PetscCall(PetscFree(dummy));
3128: PetscCall(VecRestoreArray(x, NULL));
3129: *a = NULL;
3130: PetscFunctionReturn(PETSC_SUCCESS);
3131: }
3133: /*@
3134: VecRestoreArray3dWrite - Restores a vector after `VecGetArray3dWrite()` has been called.
3136: Logically Collective
3138: Input Parameters:
3139: + x - the vector
3140: . m - first dimension of three dimensional array
3141: . n - second dimension of the three dimensional array
3142: . p - third dimension of the three dimensional array
3143: . mstart - first index you will use in first coordinate direction (often 0)
3144: . nstart - first index in the second coordinate direction (often 0)
3145: . pstart - first index in the third coordinate direction (often 0)
3146: - a - location of pointer to array obtained from VecGetArray3d()
3148: Level: developer
3150: Notes:
3151: For regular PETSc vectors this routine does not involve any copies. For
3152: any special vectors that do not store local vector data in a contiguous
3153: array, this routine will copy the data back into the underlying
3154: vector data structure from the array obtained with `VecGetArray()`.
3156: This routine actually zeros out the `a` pointer.
3158: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3159: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3160: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3161: @*/
3162: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3163: {
3164: void *dummy;
3166: PetscFunctionBegin;
3168: PetscAssertPointer(a, 8);
3170: dummy = (void *)(*a + mstart);
3171: PetscCall(PetscFree(dummy));
3172: PetscCall(VecRestoreArrayWrite(x, NULL));
3173: *a = NULL;
3174: PetscFunctionReturn(PETSC_SUCCESS);
3175: }
3177: /*@
3178: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3179: processor's portion of the vector data. You MUST call `VecRestoreArray4d()`
3180: when you no longer need access to the array.
3182: Logically Collective
3184: Input Parameters:
3185: + x - the vector
3186: . m - first dimension of four dimensional array
3187: . n - second dimension of four dimensional array
3188: . p - third dimension of four dimensional array
3189: . q - fourth dimension of four dimensional array
3190: . mstart - first index you will use in first coordinate direction (often 0)
3191: . nstart - first index in the second coordinate direction (often 0)
3192: . pstart - first index in the third coordinate direction (often 0)
3193: - qstart - first index in the fourth coordinate direction (often 0)
3195: Output Parameter:
3196: . a - location to put pointer to the array
3198: Level: developer
3200: Notes:
3201: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3202: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3203: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3204: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3206: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3208: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3209: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3210: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3211: @*/
3212: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3213: {
3214: PetscInt i, N, j, k;
3215: PetscScalar *aa, ***b, **c;
3217: PetscFunctionBegin;
3219: PetscAssertPointer(a, 10);
3221: PetscCall(VecGetLocalSize(x, &N));
3222: 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);
3223: PetscCall(VecGetArray(x, &aa));
3225: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3226: b = (PetscScalar ***)((*a) + m);
3227: c = (PetscScalar **)(b + m * n);
3228: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3229: for (i = 0; i < m; i++)
3230: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3231: for (i = 0; i < m; i++)
3232: for (j = 0; j < n; j++)
3233: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3234: *a -= mstart;
3235: PetscFunctionReturn(PETSC_SUCCESS);
3236: }
3238: /*@
3239: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3240: processor's portion of the vector data. You MUST call `VecRestoreArray4dWrite()`
3241: when you no longer need access to the array.
3243: Logically Collective
3245: Input Parameters:
3246: + x - the vector
3247: . m - first dimension of four dimensional array
3248: . n - second dimension of four dimensional array
3249: . p - third dimension of four dimensional array
3250: . q - fourth dimension of four dimensional array
3251: . mstart - first index you will use in first coordinate direction (often 0)
3252: . nstart - first index in the second coordinate direction (often 0)
3253: . pstart - first index in the third coordinate direction (often 0)
3254: - qstart - first index in the fourth coordinate direction (often 0)
3256: Output Parameter:
3257: . a - location to put pointer to the array
3259: Level: developer
3261: Notes:
3262: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3263: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3264: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3265: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3267: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3269: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3270: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3271: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3272: @*/
3273: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3274: {
3275: PetscInt i, N, j, k;
3276: PetscScalar *aa, ***b, **c;
3278: PetscFunctionBegin;
3280: PetscAssertPointer(a, 10);
3282: PetscCall(VecGetLocalSize(x, &N));
3283: 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);
3284: PetscCall(VecGetArrayWrite(x, &aa));
3286: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3287: b = (PetscScalar ***)((*a) + m);
3288: c = (PetscScalar **)(b + m * n);
3289: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3290: for (i = 0; i < m; i++)
3291: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3292: for (i = 0; i < m; i++)
3293: for (j = 0; j < n; j++)
3294: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3295: *a -= mstart;
3296: PetscFunctionReturn(PETSC_SUCCESS);
3297: }
3299: /*@
3300: VecRestoreArray4d - Restores a vector after `VecGetArray4d()` has been called.
3302: Logically Collective
3304: Input Parameters:
3305: + x - the vector
3306: . m - first dimension of four dimensional array
3307: . n - second dimension of the four dimensional array
3308: . p - third dimension of the four dimensional array
3309: . q - fourth dimension of the four dimensional array
3310: . mstart - first index you will use in first coordinate direction (often 0)
3311: . nstart - first index in the second coordinate direction (often 0)
3312: . pstart - first index in the third coordinate direction (often 0)
3313: . qstart - first index in the fourth coordinate direction (often 0)
3314: - a - location of pointer to array obtained from VecGetArray4d()
3316: Level: developer
3318: Notes:
3319: For regular PETSc vectors this routine does not involve any copies. For
3320: any special vectors that do not store local vector data in a contiguous
3321: array, this routine will copy the data back into the underlying
3322: vector data structure from the array obtained with `VecGetArray()`.
3324: This routine actually zeros out the `a` pointer.
3326: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3327: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3328: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`
3329: @*/
3330: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3331: {
3332: void *dummy;
3334: PetscFunctionBegin;
3336: PetscAssertPointer(a, 10);
3338: dummy = (void *)(*a + mstart);
3339: PetscCall(PetscFree(dummy));
3340: PetscCall(VecRestoreArray(x, NULL));
3341: *a = NULL;
3342: PetscFunctionReturn(PETSC_SUCCESS);
3343: }
3345: /*@
3346: VecRestoreArray4dWrite - Restores a vector after `VecGetArray4dWrite()` has been called.
3348: Logically Collective
3350: Input Parameters:
3351: + x - the vector
3352: . m - first dimension of four dimensional array
3353: . n - second dimension of the four dimensional array
3354: . p - third dimension of the four dimensional array
3355: . q - fourth dimension of the four dimensional array
3356: . mstart - first index you will use in first coordinate direction (often 0)
3357: . nstart - first index in the second coordinate direction (often 0)
3358: . pstart - first index in the third coordinate direction (often 0)
3359: . qstart - first index in the fourth coordinate direction (often 0)
3360: - a - location of pointer to array obtained from `VecGetArray4d()`
3362: Level: developer
3364: Notes:
3365: For regular PETSc vectors this routine does not involve any copies. For
3366: any special vectors that do not store local vector data in a contiguous
3367: array, this routine will copy the data back into the underlying
3368: vector data structure from the array obtained with `VecGetArray()`.
3370: This routine actually zeros out the `a` pointer.
3372: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3373: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3374: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3375: @*/
3376: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3377: {
3378: void *dummy;
3380: PetscFunctionBegin;
3382: PetscAssertPointer(a, 10);
3384: dummy = (void *)(*a + mstart);
3385: PetscCall(PetscFree(dummy));
3386: PetscCall(VecRestoreArrayWrite(x, NULL));
3387: *a = NULL;
3388: PetscFunctionReturn(PETSC_SUCCESS);
3389: }
3391: /*@
3392: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3393: processor's portion of the vector data. You MUST call `VecRestoreArray2dRead()`
3394: when you no longer need access to the array.
3396: Logically Collective
3398: Input Parameters:
3399: + x - the vector
3400: . m - first dimension of two dimensional array
3401: . n - second dimension of two dimensional array
3402: . mstart - first index you will use in first coordinate direction (often 0)
3403: - nstart - first index in the second coordinate direction (often 0)
3405: Output Parameter:
3406: . a - location to put pointer to the array
3408: Level: developer
3410: Notes:
3411: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
3412: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3413: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3414: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
3416: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3418: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3419: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3420: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3421: @*/
3422: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3423: {
3424: PetscInt i, N;
3425: const PetscScalar *aa;
3427: PetscFunctionBegin;
3429: PetscAssertPointer(a, 6);
3431: PetscCall(VecGetLocalSize(x, &N));
3432: 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);
3433: PetscCall(VecGetArrayRead(x, &aa));
3435: PetscCall(PetscMalloc1(m, a));
3436: for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3437: *a -= mstart;
3438: PetscFunctionReturn(PETSC_SUCCESS);
3439: }
3441: /*@
3442: VecRestoreArray2dRead - Restores a vector after `VecGetArray2dRead()` has been called.
3444: Logically Collective
3446: Input Parameters:
3447: + x - the vector
3448: . m - first dimension of two dimensional array
3449: . n - second dimension of the two dimensional array
3450: . mstart - first index you will use in first coordinate direction (often 0)
3451: . nstart - first index in the second coordinate direction (often 0)
3452: - a - location of pointer to array obtained from VecGetArray2d()
3454: Level: developer
3456: Notes:
3457: For regular PETSc vectors this routine does not involve any copies. For
3458: any special vectors that do not store local vector data in a contiguous
3459: array, this routine will copy the data back into the underlying
3460: vector data structure from the array obtained with `VecGetArray()`.
3462: This routine actually zeros out the `a` pointer.
3464: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3465: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3466: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3467: @*/
3468: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3469: {
3470: void *dummy;
3472: PetscFunctionBegin;
3474: PetscAssertPointer(a, 6);
3476: dummy = (void *)(*a + mstart);
3477: PetscCall(PetscFree(dummy));
3478: PetscCall(VecRestoreArrayRead(x, NULL));
3479: *a = NULL;
3480: PetscFunctionReturn(PETSC_SUCCESS);
3481: }
3483: /*@
3484: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3485: processor's portion of the vector data. You MUST call `VecRestoreArray1dRead()`
3486: when you no longer need access to the array.
3488: Logically Collective
3490: Input Parameters:
3491: + x - the vector
3492: . m - first dimension of two dimensional array
3493: - mstart - first index you will use in first coordinate direction (often 0)
3495: Output Parameter:
3496: . a - location to put pointer to the array
3498: Level: developer
3500: Notes:
3501: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3502: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3503: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3505: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3507: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3508: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3509: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3510: @*/
3511: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3512: {
3513: PetscInt N;
3515: PetscFunctionBegin;
3517: PetscAssertPointer(a, 4);
3519: PetscCall(VecGetLocalSize(x, &N));
3520: 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);
3521: PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3522: *a -= mstart;
3523: PetscFunctionReturn(PETSC_SUCCESS);
3524: }
3526: /*@
3527: VecRestoreArray1dRead - Restores a vector after `VecGetArray1dRead()` has been called.
3529: Logically Collective
3531: Input Parameters:
3532: + x - the vector
3533: . m - first dimension of two dimensional array
3534: . mstart - first index you will use in first coordinate direction (often 0)
3535: - a - location of pointer to array obtained from `VecGetArray1dRead()`
3537: Level: developer
3539: Notes:
3540: For regular PETSc vectors this routine does not involve any copies. For
3541: any special vectors that do not store local vector data in a contiguous
3542: array, this routine will copy the data back into the underlying
3543: vector data structure from the array obtained with `VecGetArray1dRead()`.
3545: This routine actually zeros out the `a` pointer.
3547: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3548: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3549: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3550: @*/
3551: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3552: {
3553: PetscFunctionBegin;
3556: PetscCall(VecRestoreArrayRead(x, NULL));
3557: *a = NULL;
3558: PetscFunctionReturn(PETSC_SUCCESS);
3559: }
3561: /*@
3562: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3563: processor's portion of the vector data. You MUST call `VecRestoreArray3dRead()`
3564: when you no longer need access to the array.
3566: Logically Collective
3568: Input Parameters:
3569: + x - the vector
3570: . m - first dimension of three dimensional array
3571: . n - second dimension of three dimensional array
3572: . p - third dimension of three dimensional array
3573: . mstart - first index you will use in first coordinate direction (often 0)
3574: . nstart - first index in the second coordinate direction (often 0)
3575: - pstart - first index in the third coordinate direction (often 0)
3577: Output Parameter:
3578: . a - location to put pointer to the array
3580: Level: developer
3582: Notes:
3583: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3584: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3585: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3586: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3dRead()`.
3588: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3590: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3591: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3592: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3593: @*/
3594: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3595: {
3596: PetscInt i, N, j;
3597: const PetscScalar *aa;
3598: PetscScalar **b;
3600: PetscFunctionBegin;
3602: PetscAssertPointer(a, 8);
3604: PetscCall(VecGetLocalSize(x, &N));
3605: 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);
3606: PetscCall(VecGetArrayRead(x, &aa));
3608: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3609: b = (PetscScalar **)((*a) + m);
3610: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3611: for (i = 0; i < m; i++)
3612: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset((PetscScalar *)aa, i * n * p + j * p - pstart);
3613: *a -= mstart;
3614: PetscFunctionReturn(PETSC_SUCCESS);
3615: }
3617: /*@
3618: VecRestoreArray3dRead - Restores a vector after `VecGetArray3dRead()` has been called.
3620: Logically Collective
3622: Input Parameters:
3623: + x - the vector
3624: . m - first dimension of three dimensional array
3625: . n - second dimension of the three dimensional array
3626: . p - third dimension of the three dimensional array
3627: . mstart - first index you will use in first coordinate direction (often 0)
3628: . nstart - first index in the second coordinate direction (often 0)
3629: . pstart - first index in the third coordinate direction (often 0)
3630: - a - location of pointer to array obtained from `VecGetArray3dRead()`
3632: Level: developer
3634: Notes:
3635: For regular PETSc vectors this routine does not involve any copies. For
3636: any special vectors that do not store local vector data in a contiguous
3637: array, this routine will copy the data back into the underlying
3638: vector data structure from the array obtained with `VecGetArray()`.
3640: This routine actually zeros out the `a` pointer.
3642: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3643: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3644: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3645: @*/
3646: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3647: {
3648: void *dummy;
3650: PetscFunctionBegin;
3652: PetscAssertPointer(a, 8);
3654: dummy = (void *)(*a + mstart);
3655: PetscCall(PetscFree(dummy));
3656: PetscCall(VecRestoreArrayRead(x, NULL));
3657: *a = NULL;
3658: PetscFunctionReturn(PETSC_SUCCESS);
3659: }
3661: /*@
3662: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3663: processor's portion of the vector data. You MUST call `VecRestoreArray4dRead()`
3664: when you no longer need access to the array.
3666: Logically Collective
3668: Input Parameters:
3669: + x - the vector
3670: . m - first dimension of four dimensional array
3671: . n - second dimension of four dimensional array
3672: . p - third dimension of four dimensional array
3673: . q - fourth dimension of four dimensional array
3674: . mstart - first index you will use in first coordinate direction (often 0)
3675: . nstart - first index in the second coordinate direction (often 0)
3676: . pstart - first index in the third coordinate direction (often 0)
3677: - qstart - first index in the fourth coordinate direction (often 0)
3679: Output Parameter:
3680: . a - location to put pointer to the array
3682: Level: beginner
3684: Notes:
3685: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3686: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3687: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3688: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3690: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3692: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3693: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3694: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3695: @*/
3696: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3697: {
3698: PetscInt i, N, j, k;
3699: const PetscScalar *aa;
3700: PetscScalar ***b, **c;
3702: PetscFunctionBegin;
3704: PetscAssertPointer(a, 10);
3706: PetscCall(VecGetLocalSize(x, &N));
3707: 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);
3708: PetscCall(VecGetArrayRead(x, &aa));
3710: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3711: b = (PetscScalar ***)((*a) + m);
3712: c = (PetscScalar **)(b + m * n);
3713: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3714: for (i = 0; i < m; i++)
3715: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3716: for (i = 0; i < m; i++)
3717: for (j = 0; j < n; j++)
3718: 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;
3719: *a -= mstart;
3720: PetscFunctionReturn(PETSC_SUCCESS);
3721: }
3723: /*@
3724: VecRestoreArray4dRead - Restores a vector after `VecGetArray4d()` has been called.
3726: Logically Collective
3728: Input Parameters:
3729: + x - the vector
3730: . m - first dimension of four dimensional array
3731: . n - second dimension of the four dimensional array
3732: . p - third dimension of the four dimensional array
3733: . q - fourth dimension of the four dimensional array
3734: . mstart - first index you will use in first coordinate direction (often 0)
3735: . nstart - first index in the second coordinate direction (often 0)
3736: . pstart - first index in the third coordinate direction (often 0)
3737: . qstart - first index in the fourth coordinate direction (often 0)
3738: - a - location of pointer to array obtained from `VecGetArray4dRead()`
3740: Level: beginner
3742: Notes:
3743: For regular PETSc vectors this routine does not involve any copies. For
3744: any special vectors that do not store local vector data in a contiguous
3745: array, this routine will copy the data back into the underlying
3746: vector data structure from the array obtained with `VecGetArray()`.
3748: This routine actually zeros out the `a` pointer.
3750: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3751: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3752: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3753: @*/
3754: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3755: {
3756: void *dummy;
3758: PetscFunctionBegin;
3760: PetscAssertPointer(a, 10);
3762: dummy = (void *)(*a + mstart);
3763: PetscCall(PetscFree(dummy));
3764: PetscCall(VecRestoreArrayRead(x, NULL));
3765: *a = NULL;
3766: PetscFunctionReturn(PETSC_SUCCESS);
3767: }
3769: /*@
3770: VecLockGet - Get the current lock status of a vector
3772: Logically Collective
3774: Input Parameter:
3775: . x - the vector
3777: Output Parameter:
3778: . state - greater than zero indicates the vector is locked for read; less than zero indicates the vector is
3779: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3781: Level: advanced
3783: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3784: @*/
3785: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3786: {
3787: PetscFunctionBegin;
3789: PetscAssertPointer(state, 2);
3790: *state = x->lock;
3791: PetscFunctionReturn(PETSC_SUCCESS);
3792: }
3794: /*@
3795: VecLockGetLocation - Return the source code location where a `Vec` was most recently read-locked
3797: Not Collective
3799: Input Parameter:
3800: . x - the vector
3802: Output Parameters:
3803: + file - the source file name of the most recent `VecLockReadPush()`, or `NULL` if none is active
3804: . func - the function name of the most recent `VecLockReadPush()`, or `NULL` if none is active
3805: - line - the source line number of the most recent `VecLockReadPush()`, or 0 if none is active
3807: Level: developer
3809: Note:
3810: Only produces meaningful output when PETSc is built with `PETSC_USE_DEBUG` and without threadsafety; otherwise `NULL`
3811: and 0 are returned. Intended to help debug read-lock violations reported by `VecGetArray()` and similar routines.
3813: .seealso: `Vec`, `VecLockGet()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecGetArray()`
3814: @*/
3815: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3816: {
3817: PetscFunctionBegin;
3819: PetscAssertPointer(file, 2);
3820: PetscAssertPointer(func, 3);
3821: PetscAssertPointer(line, 4);
3822: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3823: {
3824: const int index = x->lockstack.currentsize - 1;
3826: *file = index < 0 ? NULL : x->lockstack.file[index];
3827: *func = index < 0 ? NULL : x->lockstack.function[index];
3828: *line = index < 0 ? 0 : x->lockstack.line[index];
3829: }
3830: #else
3831: *file = NULL;
3832: *func = NULL;
3833: *line = 0;
3834: #endif
3835: PetscFunctionReturn(PETSC_SUCCESS);
3836: }
3838: /*@
3839: VecLockReadPush - Push a read-only lock on a vector to prevent it from being written to
3841: Logically Collective
3843: Input Parameter:
3844: . x - the vector
3846: Level: intermediate
3848: Notes:
3849: If this is set then calls to `VecGetArray()` or `VecSetValues()` or any other routines that change the vectors values will generate an error.
3851: The call can be nested, i.e., called multiple times on the same vector, but each `VecLockReadPush()` has to have one matching
3852: `VecLockReadPop()`, which removes the latest read-only lock.
3854: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3855: @*/
3856: PetscErrorCode VecLockReadPush(Vec x)
3857: {
3858: PetscFunctionBegin;
3860: 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");
3861: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3862: {
3863: const char *file, *func;
3864: int index, line;
3866: if ((index = petscstack.currentsize - 2) < 0) {
3867: // vec was locked "outside" of petsc, either in user-land or main. the error message will
3868: // now show this function as the culprit, but it will include the stacktrace
3869: file = "unknown user-file";
3870: func = "unknown_user_function";
3871: line = 0;
3872: } else {
3873: file = petscstack.file[index];
3874: func = petscstack.function[index];
3875: line = petscstack.line[index];
3876: }
3877: PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
3878: }
3879: #endif
3880: PetscFunctionReturn(PETSC_SUCCESS);
3881: }
3883: /*@
3884: VecLockReadPop - Pop a read-only lock from a vector
3886: Logically Collective
3888: Input Parameter:
3889: . x - the vector
3891: Level: intermediate
3893: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
3894: @*/
3895: PetscErrorCode VecLockReadPop(Vec x)
3896: {
3897: PetscFunctionBegin;
3899: PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
3900: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3901: {
3902: const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];
3904: PetscStackPop_Private(x->lockstack, previous);
3905: }
3906: #endif
3907: PetscFunctionReturn(PETSC_SUCCESS);
3908: }
3910: /*@
3911: VecLockWriteSet - Lock or unlock a vector for exclusive read/write access
3913: Logically Collective
3915: Input Parameters:
3916: + x - the vector
3917: - flg - `PETSC_TRUE` to lock the vector for exclusive read/write access; `PETSC_FALSE` to unlock it.
3919: Level: intermediate
3921: Notes:
3922: The function is useful in split-phase computations, which usually have a begin phase and an end phase.
3923: One can call `VecLockWriteSet`(x,`PETSC_TRUE`) in the begin phase to lock a vector for exclusive
3924: access, and call `VecLockWriteSet`(x,`PETSC_FALSE`) in the end phase to unlock the vector from exclusive
3925: access. In this way, one is ensured no other operations can access the vector in between. The code may like
3927: .vb
3928: VecGetArray(x,&xdata); // begin phase
3929: VecLockWriteSet(v,PETSC_TRUE);
3931: Other operations, which can not access x anymore (they can access xdata, of course)
3933: VecRestoreArray(x,&vdata); // end phase
3934: VecLockWriteSet(v,PETSC_FALSE);
3935: .ve
3937: The call can not be nested on the same vector, in other words, one can not call `VecLockWriteSet`(x,`PETSC_TRUE`)
3938: again before calling `VecLockWriteSet`(v,`PETSC_FALSE`).
3940: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
3941: @*/
3942: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
3943: {
3944: PetscFunctionBegin;
3946: if (flg) {
3947: 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");
3948: 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");
3949: x->lock = -1;
3950: } else {
3951: 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");
3952: x->lock = 0;
3953: }
3954: PetscFunctionReturn(PETSC_SUCCESS);
3955: }