Actual source code: rvector.c
1: /*
2: Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
3: These are the vector functions the user calls.
4: */
5: #include <petsc/private/vecimpl.h>
7: PetscInt VecGetSubVectorSavedStateId = -1;
9: #if PetscDefined(USE_DEBUG)
10: // this is a no-op '0' macro in optimized builds
11: PetscErrorCode VecValidValues_Internal(Vec vec, PetscInt argnum, PetscBool begin)
12: {
13: PetscFunctionBegin;
14: if (vec->petscnative || vec->ops->getarray) {
15: PetscInt n;
16: const PetscScalar *x;
17: PetscOffloadMask mask;
19: PetscCall(VecGetOffloadMask(vec, &mask));
20: if (!PetscOffloadHost(mask)) PetscFunctionReturn(PETSC_SUCCESS);
21: PetscCall(VecGetLocalSize(vec, &n));
22: PetscCall(VecGetArrayRead(vec, &x));
23: for (PetscInt i = 0; i < n; i++) {
24: if (begin) {
25: PetscCheck(!PetscIsInfOrNanScalar(x[i]), PETSC_COMM_SELF, PETSC_ERR_FP, "Vec entry at local location %" PetscInt_FMT " is not-a-number or infinite at beginning of function: Parameter number %" PetscInt_FMT, i, argnum);
26: } else {
27: PetscCheck(!PetscIsInfOrNanScalar(x[i]), PETSC_COMM_SELF, PETSC_ERR_FP, "Vec entry at local location %" PetscInt_FMT " is not-a-number or infinite at end of function: Parameter number %" PetscInt_FMT, i, argnum);
28: }
29: }
30: PetscCall(VecRestoreArrayRead(vec, &x));
31: }
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
34: #endif
36: /*@
37: VecMaxPointwiseDivide - Computes the maximum of the componentwise division `max = max_i abs(x[i]/y[i])`.
39: Logically Collective
41: Input Parameters:
42: + x - the numerators
43: - y - the denominators
45: Output Parameter:
46: . max - the result
48: Level: advanced
50: Notes:
51: `x` and `y` may be the same vector
53: if a particular `y[i]` is zero, it is treated as 1 in the above formula
55: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`
56: @*/
57: PetscErrorCode VecMaxPointwiseDivide(Vec x, Vec y, PetscReal *max)
58: {
59: PetscFunctionBegin;
62: PetscAssertPointer(max, 3);
65: PetscCheckSameTypeAndComm(x, 1, y, 2);
66: VecCheckSameSize(x, 1, y, 2);
67: VecCheckAssembled(x);
68: VecCheckAssembled(y);
69: PetscCall(VecLockReadPush(x));
70: PetscCall(VecLockReadPush(y));
71: PetscUseTypeMethod(x, maxpointwisedivide, y, max);
72: PetscCall(VecLockReadPop(x));
73: PetscCall(VecLockReadPop(y));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: /*@
78: VecDot - Computes the vector dot product.
80: Collective
82: Input Parameters:
83: + x - first vector
84: - y - second vector
86: Output Parameter:
87: . val - the dot product
89: Level: intermediate
91: Notes for Users of Complex Numbers:
92: For complex vectors, `VecDot()` computes
93: .vb
94: val = (x,y) = y^H x,
95: .ve
96: where y^H denotes the conjugate transpose of y. Note that this corresponds to the usual "mathematicians" complex
97: inner product where the SECOND argument gets the complex conjugate. Since the `BLASdot()` complex conjugates the first
98: first argument we call the `BLASdot()` with the arguments reversed.
100: Use `VecTDot()` for the indefinite form
101: .vb
102: val = (x,y) = y^T x,
103: .ve
104: where y^T denotes the transpose of y.
106: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDotRealPart()`
107: @*/
108: PetscErrorCode VecDot(Vec x, Vec y, PetscScalar *val)
109: {
110: PetscFunctionBegin;
113: PetscAssertPointer(val, 3);
116: PetscCheckSameTypeAndComm(x, 1, y, 2);
117: VecCheckSameSize(x, 1, y, 2);
118: VecCheckAssembled(x);
119: VecCheckAssembled(y);
121: PetscCall(VecLockReadPush(x));
122: PetscCall(VecLockReadPush(y));
123: PetscCall(PetscLogEventBegin(VEC_Dot, x, y, 0, 0));
124: PetscUseTypeMethod(x, dot, y, val);
125: PetscCall(PetscLogEventEnd(VEC_Dot, x, y, 0, 0));
126: PetscCall(VecLockReadPop(x));
127: PetscCall(VecLockReadPop(y));
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: /*@
132: VecDotRealPart - Computes the real part of the vector dot product.
134: Collective
136: Input Parameters:
137: + x - first vector
138: - y - second vector
140: Output Parameter:
141: . val - the real part of the dot product;
143: Level: intermediate
145: Notes for Users of Complex Numbers:
146: See `VecDot()` for more details on the definition of the dot product for complex numbers
148: For real numbers this returns the same value as `VecDot()`
150: For complex numbers in C^n (that is a vector of n components with a complex number for each component) this is equal to the usual real dot product on the
151: the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)
153: Developer Notes:
154: This is not currently optimized to compute only the real part of the dot product.
156: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDot()`, `VecDotNorm2()`
157: @*/
158: PetscErrorCode VecDotRealPart(Vec x, Vec y, PetscReal *val)
159: {
160: PetscScalar fdot;
162: PetscFunctionBegin;
163: PetscCall(VecDot(x, y, &fdot));
164: *val = PetscRealPart(fdot);
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: /*@
169: VecNorm - Computes the vector norm.
171: Collective
173: Input Parameters:
174: + x - the vector
175: - type - the type of the norm requested
177: Output Parameter:
178: . val - the norm
180: Level: intermediate
182: Notes:
183: See `NormType` for descriptions of each norm.
185: For complex numbers `NORM_1` will return the traditional 1 norm of the 2 norm of the complex
186: numbers; that is the 1 norm of the absolute values of the complex entries. In PETSc 3.6 and
187: earlier releases it returned the 1 norm of the 1 norm of the complex entries (what is
188: returned by the BLAS routine `asum()`). Both are valid norms but most people expect the former.
190: This routine stashes the computed norm value, repeated calls before the vector entries are
191: changed are then rapid since the precomputed value is immediately available. Certain vector
192: operations such as `VecSet()` store the norms so the value is immediately available and does
193: not need to be explicitly computed. `VecScale()` updates any stashed norm values, thus calls
194: after `VecScale()` do not need to explicitly recompute the norm.
196: .seealso: [](ch_vectors), `Vec`, `NormType`, `VecDot()`, `VecTDot()`, `VecDotBegin()`, `VecDotEnd()`, `VecNormAvailable()`,
197: `VecNormBegin()`, `VecNormEnd()`, `NormType()`
198: @*/
199: PetscErrorCode VecNorm(Vec x, NormType type, PetscReal *val)
200: {
201: PetscBool flg = PETSC_TRUE;
203: PetscFunctionBegin;
206: VecCheckAssembled(x);
208: PetscAssertPointer(val, 3);
210: PetscCall(VecNormAvailable(x, type, &flg, val));
211: // check that all MPI processes call this routine together and have same availability
212: if (PetscDefined(USE_DEBUG)) {
213: PetscMPIInt b0 = (PetscMPIInt)flg, b1[2], b2[2];
214: b1[0] = -b0;
215: b1[1] = b0;
216: PetscCallMPI(MPIU_Allreduce(b1, b2, 2, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
217: PetscCheck(-b2[0] == b2[1], PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_WRONGSTATE, "Some MPI processes have cached %s norm, others do not. This may happen when some MPI processes call VecGetArray() and some others do not.", NormTypes[type]);
218: if (flg) {
219: PetscReal b1[2], b2[2];
220: b1[0] = -(*val);
221: b1[1] = *val;
222: PetscCallMPI(MPIU_Allreduce(b1, b2, 2, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)x)));
223: PetscCheck((PetscIsNanReal(b2[0]) && PetscIsNanReal(b2[1])) || (-b2[0] == b2[1]), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Difference in cached %s norms: local %g", NormTypes[type], (double)*val);
224: }
225: }
226: if (flg) PetscFunctionReturn(PETSC_SUCCESS);
228: PetscCall(VecLockReadPush(x));
229: PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
230: PetscUseTypeMethod(x, norm, type, val);
231: PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
232: PetscCall(VecLockReadPop(x));
234: if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: /*@
239: VecNormAvailable - Returns the vector norm if it is already known. That is, it has been previously computed and cached in the vector
241: Not Collective
243: Input Parameters:
244: + x - the vector
245: - type - one of `NORM_1` (sum_i |x[i]|), `NORM_2` sqrt(sum_i (x[i])^2), `NORM_INFINITY` max_i |x[i]|. Also available
246: `NORM_1_AND_2`, which computes both norms and stores them
247: in a two element array.
249: Output Parameters:
250: + available - `PETSC_TRUE` if the val returned is valid
251: - val - the norm
253: Level: intermediate
255: Developer Notes:
256: `PETSC_HAVE_SLOW_BLAS_NORM2` will cause a C (loop unrolled) version of the norm to be used, rather
257: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
258: (as opposed to vendor provided) because the FORTRAN BLAS `NRM2()` routine is very slow.
260: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
261: `VecNormBegin()`, `VecNormEnd()`
262: @*/
263: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
264: {
265: PetscFunctionBegin;
268: PetscAssertPointer(available, 3);
269: PetscAssertPointer(val, 4);
271: if (type == NORM_1_AND_2) {
272: *available = PETSC_FALSE;
273: } else {
274: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
275: }
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /*@
280: VecNormalize - Normalizes a vector by its 2-norm.
282: Collective
284: Input Parameter:
285: . x - the vector
287: Output Parameter:
288: . val - the vector norm before normalization. May be `NULL` if the value is not needed.
290: Level: intermediate
292: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
293: @*/
294: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
295: {
296: PetscReal norm;
298: PetscFunctionBegin;
301: PetscCall(VecSetErrorIfLocked(x, 1));
302: if (val) PetscAssertPointer(val, 2);
303: PetscCall(PetscLogEventBegin(VEC_Normalize, x, 0, 0, 0));
304: PetscCall(VecNorm(x, NORM_2, &norm));
305: if (norm == 0.0) PetscCall(PetscInfo(x, "Vector of zero norm can not be normalized; Returning only the zero norm\n"));
306: else if (PetscIsInfOrNanReal(norm)) PetscCall(PetscInfo(x, "Vector with Inf or Nan norm can not be normalized; Returning only the norm\n"));
307: else {
308: PetscScalar s = 1.0 / norm;
309: PetscCall(VecScale(x, s));
310: }
311: PetscCall(PetscLogEventEnd(VEC_Normalize, x, 0, 0, 0));
312: if (val) *val = norm;
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: /*@
317: VecMax - Determines the vector component with maximum real part and its location.
319: Collective
321: Input Parameter:
322: . x - the vector
324: Output Parameters:
325: + p - the index of `val` (pass `NULL` if you don't want this) in the vector
326: - val - the maximum component
328: Level: intermediate
330: Notes:
331: Returns the value `PETSC_MIN_REAL` and negative `p` if the vector is of length 0.
333: Returns the smallest index with the maximum value
335: Developer Note:
336: The Nag Fortran compiler does not like the symbol name VecMax
338: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
339: @*/
340: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
341: {
342: PetscFunctionBegin;
345: VecCheckAssembled(x);
346: if (p) PetscAssertPointer(p, 2);
347: PetscAssertPointer(val, 3);
348: PetscCall(VecLockReadPush(x));
349: PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
350: PetscUseTypeMethod(x, max, p, val);
351: PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
352: PetscCall(VecLockReadPop(x));
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: /*@
357: VecMin - Determines the vector component with minimum real part and its location.
359: Collective
361: Input Parameter:
362: . x - the vector
364: Output Parameters:
365: + p - the index of `val` (pass `NULL` if you don't want this location) in the vector
366: - val - the minimum component
368: Level: intermediate
370: Notes:
371: Returns the value `PETSC_MAX_REAL` and negative `p` if the vector is of length 0.
373: This returns the smallest index with the minimum value
375: Developer Note:
376: The Nag Fortran compiler does not like the symbol name VecMin
378: .seealso: [](ch_vectors), `Vec`, `VecMax()`
379: @*/
380: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
381: {
382: PetscFunctionBegin;
385: VecCheckAssembled(x);
386: if (p) PetscAssertPointer(p, 2);
387: PetscAssertPointer(val, 3);
388: PetscCall(VecLockReadPush(x));
389: PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
390: PetscUseTypeMethod(x, min, p, val);
391: PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
392: PetscCall(VecLockReadPop(x));
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: /*@
397: VecTDot - Computes an indefinite vector dot product. That is, this
398: routine does NOT use the complex conjugate.
400: Collective
402: Input Parameters:
403: + x - first vector
404: - y - second vector
406: Output Parameter:
407: . val - the dot product
409: Level: intermediate
411: Notes for Users of Complex Numbers:
412: For complex vectors, `VecTDot()` computes the indefinite form
413: .vb
414: val = (x,y) = y^T x,
415: .ve
416: where y^T denotes the transpose of y.
418: Use `VecDot()` for the inner product
419: .vb
420: val = (x,y) = y^H x,
421: .ve
422: where y^H denotes the conjugate transpose of y.
424: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
425: @*/
426: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
427: {
428: PetscFunctionBegin;
431: PetscAssertPointer(val, 3);
434: PetscCheckSameTypeAndComm(x, 1, y, 2);
435: VecCheckSameSize(x, 1, y, 2);
436: VecCheckAssembled(x);
437: VecCheckAssembled(y);
439: PetscCall(VecLockReadPush(x));
440: PetscCall(VecLockReadPush(y));
441: PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
442: PetscUseTypeMethod(x, tdot, y, val);
443: PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
444: PetscCall(VecLockReadPop(x));
445: PetscCall(VecLockReadPop(y));
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: PetscErrorCode VecScaleAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
450: {
451: PetscReal norms[4];
452: PetscBool flgs[4];
453: PetscScalar one = 1.0;
455: PetscFunctionBegin;
458: VecCheckAssembled(x);
459: PetscCall(VecSetErrorIfLocked(x, 1));
461: if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);
463: /* get current stashed norms */
464: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
466: PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
467: VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
468: PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));
470: PetscCall(PetscObjectStateIncrease((PetscObject)x));
471: /* put the scaled stashed norms back into the Vec */
472: for (PetscInt i = 0; i < 4; i++) {
473: PetscReal ar = PetscAbsScalar(alpha);
474: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
475: }
476: PetscFunctionReturn(PETSC_SUCCESS);
477: }
479: /*@
480: VecScale - Scales a vector.
482: Logically Collective
484: Input Parameters:
485: + x - the vector
486: - alpha - the scalar
488: Level: intermediate
490: Note:
491: For a vector with n components, `VecScale()` computes x[i] = alpha * x[i], for i=1,...,n.
493: .seealso: [](ch_vectors), `Vec`, `VecSet()`
494: @*/
495: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
496: {
497: PetscFunctionBegin;
498: PetscCall(VecScaleAsync_Private(x, alpha, NULL));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
503: {
504: PetscFunctionBegin;
507: VecCheckAssembled(x);
509: PetscCall(VecSetErrorIfLocked(x, 1));
511: if (alpha == 0) {
512: PetscReal norm;
513: PetscBool set;
515: PetscCall(VecNormAvailable(x, NORM_2, &set, &norm));
516: if (set == PETSC_TRUE && norm == 0) PetscFunctionReturn(PETSC_SUCCESS);
517: }
518: PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
519: VecMethodDispatch(x, dctx, VecAsyncFnName(Set), set, (Vec, PetscScalar, PetscDeviceContext), alpha);
520: PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
521: PetscCall(PetscObjectStateIncrease((PetscObject)x));
523: /* norms can be simply set (if |alpha|*N not too large) */
524: {
525: PetscReal val = PetscAbsScalar(alpha);
526: const PetscInt N = x->map->N;
528: if (N == 0) {
529: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0));
530: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
531: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
532: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
533: } else if (val > PETSC_MAX_REAL / N) {
534: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
535: } else {
536: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
537: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
538: val *= PetscSqrtReal((PetscReal)N);
539: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
540: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
541: }
542: }
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: /*@
547: VecSet - Sets all components of a vector to a single scalar value.
549: Logically Collective
551: Input Parameters:
552: + x - the vector
553: - alpha - the scalar
555: Level: beginner
557: Notes:
558: For a vector of dimension n, `VecSet()` sets x[i] = alpha, for i=1,...,n,
559: so that all vector entries then equal the identical
560: scalar value, `alpha`. Use the more general routine
561: `VecSetValues()` to set different vector entries.
563: You CANNOT call this after you have called `VecSetValues()` but before you call
564: `VecAssemblyBegin()`
566: If `alpha` is zero and the norm of the vector is known to be zero then this skips the unneeded zeroing process
568: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
569: @*/
570: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
571: {
572: PetscFunctionBegin;
573: PetscCall(VecSetAsync_Private(x, alpha, NULL));
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: PetscErrorCode VecAXPYAsync_Private(Vec y, PetscScalar alpha, Vec x, PetscDeviceContext dctx)
578: {
579: PetscFunctionBegin;
584: PetscCheckSameTypeAndComm(x, 3, y, 1);
585: VecCheckSameSize(x, 3, y, 1);
586: VecCheckAssembled(x);
587: VecCheckAssembled(y);
589: if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
590: PetscCall(VecSetErrorIfLocked(y, 1));
591: if (x == y) {
592: PetscCall(VecScale(y, alpha + 1.0));
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
595: PetscCall(VecLockReadPush(x));
596: PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
597: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPY), axpy, (Vec, PetscScalar, Vec, PetscDeviceContext), alpha, x);
598: PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
599: PetscCall(VecLockReadPop(x));
600: PetscCall(PetscObjectStateIncrease((PetscObject)y));
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }
603: /*@
604: VecAXPY - Computes `y = alpha x + y`.
606: Logically Collective
608: Input Parameters:
609: + alpha - the scalar
610: . x - vector scale by `alpha`
611: - y - vector accumulated into
613: Output Parameter:
614: . y - output vector
616: Level: intermediate
618: Notes:
619: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
620: .vb
621: VecAXPY(y,alpha,x) y = alpha x + y
622: VecAYPX(y,beta,x) y = x + beta y
623: VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
624: VecWAXPY(w,alpha,x,y) w = alpha x + y
625: VecAXPBYPCZ(z,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
626: VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
627: .ve
629: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
630: @*/
631: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
632: {
633: PetscFunctionBegin;
634: PetscCall(VecAXPYAsync_Private(y, alpha, x, NULL));
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: PetscErrorCode VecAYPXAsync_Private(Vec y, PetscScalar beta, Vec x, PetscDeviceContext dctx)
639: {
640: PetscFunctionBegin;
645: PetscCheckSameTypeAndComm(x, 3, y, 1);
646: VecCheckSameSize(x, 1, y, 3);
647: VecCheckAssembled(x);
648: VecCheckAssembled(y);
650: PetscCall(VecSetErrorIfLocked(y, 1));
651: if (x == y) {
652: PetscCall(VecScale(y, beta + 1.0));
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
655: PetscCall(VecLockReadPush(x));
656: if (beta == (PetscScalar)0.0) {
657: PetscCall(VecCopy(x, y));
658: } else {
659: PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
660: VecMethodDispatch(y, dctx, VecAsyncFnName(AYPX), aypx, (Vec, PetscScalar, Vec, PetscDeviceContext), beta, x);
661: PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
662: PetscCall(PetscObjectStateIncrease((PetscObject)y));
663: }
664: PetscCall(VecLockReadPop(x));
665: PetscFunctionReturn(PETSC_SUCCESS);
666: }
668: /*@
669: VecAYPX - Computes `y = x + beta y`.
671: Logically Collective
673: Input Parameters:
674: + beta - the scalar
675: . x - the unscaled vector
676: - y - the vector to be scaled
678: Output Parameter:
679: . y - output vector
681: Level: intermediate
683: Developer Notes:
684: The implementation is optimized for `beta` of -1.0, 0.0, and 1.0
686: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
687: @*/
688: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
689: {
690: PetscFunctionBegin;
691: PetscCall(VecAYPXAsync_Private(y, beta, x, NULL));
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: PetscErrorCode VecAXPBYAsync_Private(Vec y, PetscScalar alpha, PetscScalar beta, Vec x, PetscDeviceContext dctx)
696: {
697: PetscFunctionBegin;
702: PetscCheckSameTypeAndComm(x, 4, y, 1);
703: VecCheckSameSize(y, 1, x, 4);
704: VecCheckAssembled(x);
705: VecCheckAssembled(y);
708: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
709: if (x == y) {
710: PetscCall(VecScale(y, alpha + beta));
711: PetscFunctionReturn(PETSC_SUCCESS);
712: }
714: PetscCall(VecSetErrorIfLocked(y, 1));
715: PetscCall(VecLockReadPush(x));
716: PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
717: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPBY), axpby, (Vec, PetscScalar, PetscScalar, Vec, PetscDeviceContext), alpha, beta, x);
718: PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
719: PetscCall(PetscObjectStateIncrease((PetscObject)y));
720: PetscCall(VecLockReadPop(x));
721: PetscFunctionReturn(PETSC_SUCCESS);
722: }
724: /*@
725: VecAXPBY - Computes `y = alpha x + beta y`.
727: Logically Collective
729: Input Parameters:
730: + alpha - first scalar
731: . beta - second scalar
732: . x - the first scaled vector
733: - y - the second scaled vector
735: Output Parameter:
736: . y - output vector
738: Level: intermediate
740: Developer Notes:
741: The implementation is optimized for `alpha` and/or `beta` values of 0.0 and 1.0
743: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
744: @*/
745: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
746: {
747: PetscFunctionBegin;
748: PetscCall(VecAXPBYAsync_Private(y, alpha, beta, x, NULL));
749: PetscFunctionReturn(PETSC_SUCCESS);
750: }
752: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
753: {
754: PetscFunctionBegin;
761: PetscCheckSameTypeAndComm(x, 5, y, 6);
762: PetscCheckSameTypeAndComm(x, 5, z, 1);
763: VecCheckSameSize(x, 5, y, 6);
764: VecCheckSameSize(x, 5, z, 1);
765: PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
766: PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
767: VecCheckAssembled(x);
768: VecCheckAssembled(y);
769: VecCheckAssembled(z);
773: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
775: PetscCall(VecSetErrorIfLocked(z, 1));
776: PetscCall(VecLockReadPush(x));
777: PetscCall(VecLockReadPush(y));
778: PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
779: VecMethodDispatch(z, dctx, VecAsyncFnName(AXPBYPCZ), axpbypcz, (Vec, PetscScalar, PetscScalar, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, beta, gamma, x, y);
780: PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
781: PetscCall(PetscObjectStateIncrease((PetscObject)z));
782: PetscCall(VecLockReadPop(x));
783: PetscCall(VecLockReadPop(y));
784: PetscFunctionReturn(PETSC_SUCCESS);
785: }
786: /*@
787: VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`
789: Logically Collective
791: Input Parameters:
792: + alpha - first scalar
793: . beta - second scalar
794: . gamma - third scalar
795: . x - first vector
796: . y - second vector
797: - z - third vector
799: Output Parameter:
800: . z - output vector
802: Level: intermediate
804: Note:
805: `x`, `y` and `z` must be different vectors
807: Developer Notes:
808: The implementation is optimized for `alpha` of 1.0 and `gamma` of 1.0 or 0.0
810: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
811: @*/
812: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
813: {
814: PetscFunctionBegin;
815: PetscCall(VecAXPBYPCZAsync_Private(z, alpha, beta, gamma, x, y, NULL));
816: PetscFunctionReturn(PETSC_SUCCESS);
817: }
819: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
820: {
821: PetscFunctionBegin;
828: PetscCheckSameTypeAndComm(x, 3, y, 4);
829: PetscCheckSameTypeAndComm(y, 4, w, 1);
830: VecCheckSameSize(x, 3, y, 4);
831: VecCheckSameSize(x, 3, w, 1);
832: PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
833: PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
834: VecCheckAssembled(x);
835: VecCheckAssembled(y);
837: PetscCall(VecSetErrorIfLocked(w, 1));
839: PetscCall(VecLockReadPush(x));
840: PetscCall(VecLockReadPush(y));
841: if (alpha == (PetscScalar)0.0) {
842: PetscCall(VecCopyAsync_Private(y, w, dctx));
843: } else {
844: PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
845: VecMethodDispatch(w, dctx, VecAsyncFnName(WAXPY), waxpy, (Vec, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, x, y);
846: PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
847: PetscCall(PetscObjectStateIncrease((PetscObject)w));
848: }
849: PetscCall(VecLockReadPop(x));
850: PetscCall(VecLockReadPop(y));
851: PetscFunctionReturn(PETSC_SUCCESS);
852: }
854: /*@
855: VecWAXPY - Computes `w = alpha x + y`.
857: Logically Collective
859: Input Parameters:
860: + alpha - the scalar
861: . x - first vector, multiplied by `alpha`
862: - y - second vector
864: Output Parameter:
865: . w - the result
867: Level: intermediate
869: Note:
870: `w` cannot be either `x` or `y`, but `x` and `y` can be the same
872: Developer Notes:
873: The implementation is optimized for alpha of -1.0, 0.0, and 1.0
875: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
876: @*/
877: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
878: {
879: PetscFunctionBegin;
880: PetscCall(VecWAXPYAsync_Private(w, alpha, x, y, NULL));
881: PetscFunctionReturn(PETSC_SUCCESS);
882: }
884: /*@
885: VecSetValues - Inserts or adds values into certain locations of a vector.
887: Not Collective
889: Input Parameters:
890: + x - vector to insert in
891: . ni - number of elements to add
892: . ix - indices where to add
893: . y - array of values. Pass `NULL` to set all zeroes.
894: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries
896: Level: beginner
898: Notes:
899: .vb
900: `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
901: .ve
903: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
904: options cannot be mixed without intervening calls to the assembly
905: routines.
907: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
908: MUST be called after all calls to `VecSetValues()` have been completed.
910: VecSetValues() uses 0-based indices in Fortran as well as in C.
912: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
913: negative indices may be passed in ix. These rows are
914: simply ignored. This allows easily inserting element load matrices
915: with homogeneous Dirichlet boundary conditions that you don't want represented
916: in the vector.
918: Fortran Note:
919: If any of `ix` and `y` are scalars pass them using, for example,
920: .vb
921: call VecSetValues(mat, one, [ix], [y], INSERT_VALUES, ierr)
922: .ve
924: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
925: `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
926: @*/
927: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
928: {
929: PetscFunctionBeginHot;
931: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
932: PetscAssertPointer(ix, 3);
933: if (y) PetscAssertPointer(y, 4);
936: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
937: PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
938: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
939: PetscCall(PetscObjectStateIncrease((PetscObject)x));
940: PetscFunctionReturn(PETSC_SUCCESS);
941: }
943: /*@
944: VecGetValues - Gets values from certain locations of a vector. Currently
945: can only get values on the same processor on which they are owned
947: Not Collective
949: Input Parameters:
950: + x - vector to get values from
951: . ni - number of elements to get
952: - ix - indices where to get them from (in global 1d numbering)
954: Output Parameter:
955: . y - array of values, must be passed in with a length of `ni`
957: Level: beginner
959: Notes:
960: The user provides the allocated array y; it is NOT allocated in this routine
962: `VecGetValues()` gets y[i] = x[ix[i]], for i=0,...,ni-1.
964: `VecAssemblyBegin()` and `VecAssemblyEnd()` MUST be called before calling this if `VecSetValues()` or related routine has been called
966: VecGetValues() uses 0-based indices in Fortran as well as in C.
968: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
969: negative indices may be passed in ix. These rows are
970: simply ignored.
972: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
973: @*/
974: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
975: {
976: PetscFunctionBegin;
978: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
979: PetscAssertPointer(ix, 3);
980: PetscAssertPointer(y, 4);
982: VecCheckAssembled(x);
983: PetscUseTypeMethod(x, getvalues, ni, ix, y);
984: PetscFunctionReturn(PETSC_SUCCESS);
985: }
987: /*@
988: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
990: Not Collective
992: Input Parameters:
993: + x - vector to insert in
994: . ni - number of blocks to add
995: . ix - indices where to add in block count, rather than element count
996: . y - array of values. Pass `NULL` to set all zeroes.
997: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries
999: Level: intermediate
1001: Notes:
1002: `VecSetValuesBlocked()` sets x[bs*ix[i]+j] = y[bs*i+j],
1003: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
1005: Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
1006: options cannot be mixed without intervening calls to the assembly
1007: routines.
1009: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1010: MUST be called after all calls to `VecSetValuesBlocked()` have been completed.
1012: `VecSetValuesBlocked()` uses 0-based indices in Fortran as well as in C.
1014: Negative indices may be passed in ix, these rows are
1015: simply ignored. This allows easily inserting element load matrices
1016: with homogeneous Dirichlet boundary conditions that you don't want represented
1017: in the vector.
1019: Fortran Note:
1020: If any of `ix` and `y` are scalars pass them using, for example,
1021: .vb
1022: call VecSetValuesBlocked(mat, one, [ix], [y], INSERT_VALUES, ierr)
1023: .ve
1025: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
1026: `VecSetValues()`
1027: @*/
1028: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1029: {
1030: PetscFunctionBeginHot;
1032: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1033: PetscAssertPointer(ix, 3);
1034: if (y) PetscAssertPointer(y, 4);
1037: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1038: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1039: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1040: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1041: PetscFunctionReturn(PETSC_SUCCESS);
1042: }
1044: /*@
1045: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1046: using a local ordering of the nodes.
1048: Not Collective
1050: Input Parameters:
1051: + x - vector to insert in
1052: . ni - number of elements to add
1053: . ix - indices where to add
1054: . y - array of values. Pass `NULL` to set all zeroes.
1055: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1057: Level: intermediate
1059: Notes:
1060: `VecSetValuesLocal()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
1062: Calls to `VecSetValuesLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1063: options cannot be mixed without intervening calls to the assembly
1064: routines.
1066: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1067: MUST be called after all calls to `VecSetValuesLocal()` have been completed.
1069: `VecSetValuesLocal()` uses 0-based indices in Fortran as well as in C.
1071: Fortran Note:
1072: If any of `ix` and `y` are scalars pass them using, for example,
1073: .vb
1074: call VecSetValuesLocal(mat, one, [ix], [y], INSERT_VALUES, ierr)
1075: .ve
1077: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1078: `VecSetValuesBlockedLocal()`
1079: @*/
1080: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1081: {
1082: PetscInt lixp[128], *lix = lixp;
1084: PetscFunctionBeginHot;
1086: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1087: PetscAssertPointer(ix, 3);
1088: if (y) PetscAssertPointer(y, 4);
1091: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1092: if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
1093: if (x->map->mapping) {
1094: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1095: PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1096: PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1097: if (ni > 128) PetscCall(PetscFree(lix));
1098: } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1099: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1100: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1101: PetscFunctionReturn(PETSC_SUCCESS);
1102: }
1104: /*@
1105: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1106: using a local ordering of the nodes.
1108: Not Collective
1110: Input Parameters:
1111: + x - vector to insert in
1112: . ni - number of blocks to add
1113: . ix - indices where to add in block count, not element count
1114: . y - array of values. Pass `NULL` to set all zeroes.
1115: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1117: Level: intermediate
1119: Notes:
1120: `VecSetValuesBlockedLocal()` sets x[bs*ix[i]+j] = y[bs*i+j],
1121: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with `VecSetBlockSize()`.
1123: Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1124: options cannot be mixed without intervening calls to the assembly
1125: routines.
1127: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1128: MUST be called after all calls to `VecSetValuesBlockedLocal()` have been completed.
1130: `VecSetValuesBlockedLocal()` uses 0-based indices in Fortran as well as in C.
1132: Fortran Note:
1133: If any of `ix` and `y` are scalars pass them using, for example,
1134: .vb
1135: call VecSetValuesBlockedLocal(mat, one, [ix], [y], INSERT_VALUES, ierr)
1136: .ve
1138: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1139: `VecSetLocalToGlobalMapping()`
1140: @*/
1141: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1142: {
1143: PetscInt lixp[128], *lix = lixp;
1145: PetscFunctionBeginHot;
1147: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1148: PetscAssertPointer(ix, 3);
1149: if (y) PetscAssertPointer(y, 4);
1151: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1152: if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
1153: if (x->map->mapping) {
1154: if (ni > (PetscInt)PETSC_STATIC_ARRAY_LENGTH(lixp)) PetscCall(PetscMalloc1(ni, &lix));
1155: PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1156: PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1157: if (ni > (PetscInt)PETSC_STATIC_ARRAY_LENGTH(lixp)) PetscCall(PetscFree(lix));
1158: } else {
1159: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1160: }
1161: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1162: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1163: PetscFunctionReturn(PETSC_SUCCESS);
1164: }
1166: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1167: {
1168: PetscFunctionBegin;
1171: VecCheckAssembled(x);
1173: if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1174: PetscAssertPointer(y, 3);
1175: for (PetscInt i = 0; i < nv; ++i) {
1178: PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1179: VecCheckSameSize(x, 1, y[i], 3);
1180: VecCheckAssembled(y[i]);
1181: PetscCall(VecLockReadPush(y[i]));
1182: }
1183: PetscAssertPointer(result, 4);
1186: PetscCall(VecLockReadPush(x));
1187: PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1188: PetscCall((*mxdot)(x, nv, y, result));
1189: PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1190: PetscCall(VecLockReadPop(x));
1191: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1192: PetscFunctionReturn(PETSC_SUCCESS);
1193: }
1195: /*@
1196: VecMTDot - Computes indefinite vector multiple dot products.
1197: That is, it does NOT use the complex conjugate.
1199: Collective
1201: Input Parameters:
1202: + x - one vector
1203: . nv - number of vectors
1204: - y - array of vectors. Note that vectors are pointers
1206: Output Parameter:
1207: . val - array of the dot products
1209: Level: intermediate
1211: Notes for Users of Complex Numbers:
1212: For complex vectors, `VecMTDot()` computes the indefinite form
1213: .vb
1214: val = (x,y) = y^T x,
1215: .ve
1216: where y^T denotes the transpose of y.
1218: Use `VecMDot()` for the inner product
1219: .vb
1220: val = (x,y) = y^H x,
1221: .ve
1222: where y^H denotes the conjugate transpose of y.
1224: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1225: @*/
1226: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1227: {
1228: PetscFunctionBegin;
1230: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1231: PetscFunctionReturn(PETSC_SUCCESS);
1232: }
1234: /*@
1235: VecMDot - Computes multiple vector dot products.
1237: Collective
1239: Input Parameters:
1240: + x - one vector
1241: . nv - number of vectors
1242: - y - array of vectors.
1244: Output Parameter:
1245: . val - array of the dot products (does not allocate the array)
1247: Level: intermediate
1249: Notes for Users of Complex Numbers:
1250: For complex vectors, `VecMDot()` computes
1251: .vb
1252: val = (x,y) = y^H x,
1253: .ve
1254: where y^H denotes the conjugate transpose of y.
1256: Use `VecMTDot()` for the indefinite form
1257: .vb
1258: val = (x,y) = y^T x,
1259: .ve
1260: where y^T denotes the transpose of y.
1262: Note:
1263: The implementation may use BLAS 2 operations when the vectors `y` have been obtained with `VecDuplicateVecs()`
1265: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`, `VecDuplicateVecs()`
1266: @*/
1267: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1268: {
1269: PetscFunctionBegin;
1271: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1272: PetscFunctionReturn(PETSC_SUCCESS);
1273: }
1275: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
1276: {
1277: PetscFunctionBegin;
1279: VecCheckAssembled(y);
1281: PetscCall(VecSetErrorIfLocked(y, 1));
1282: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1283: if (nv) {
1284: PetscInt zeros = 0;
1286: PetscAssertPointer(alpha, 3);
1287: PetscAssertPointer(x, 4);
1288: for (PetscInt i = 0; i < nv; ++i) {
1292: PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1293: VecCheckSameSize(y, 1, x[i], 4);
1294: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1295: VecCheckAssembled(x[i]);
1296: PetscCall(VecLockReadPush(x[i]));
1297: zeros += alpha[i] == (PetscScalar)0.0;
1298: }
1300: if (zeros < nv) {
1301: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1302: VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1303: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1304: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1305: }
1307: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1308: }
1309: PetscFunctionReturn(PETSC_SUCCESS);
1310: }
1312: /*@
1313: VecMAXPY - Computes `y = y + sum alpha[i] x[i]`
1315: Logically Collective
1317: Input Parameters:
1318: + nv - number of scalars and `x` vectors
1319: . alpha - array of scalars
1320: . y - one vector
1321: - x - array of vectors
1323: Level: intermediate
1325: Notes:
1326: `y` cannot be any of the `x` vectors
1328: The implementation may use BLAS 2 operations when the vectors `y` have been obtained with `VecDuplicateVecs()`
1330: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`,`VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`, `VecDuplicateVecs()`
1331: @*/
1332: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1333: {
1334: PetscFunctionBegin;
1335: PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1336: PetscFunctionReturn(PETSC_SUCCESS);
1337: }
1339: /*@
1340: VecMAXPBY - Computes `y = beta y + sum alpha[i] x[i]`
1342: Logically Collective
1344: Input Parameters:
1345: + nv - number of scalars and `x` vectors
1346: . alpha - array of scalars
1347: . beta - scalar
1348: . y - one vector
1349: - x - array of vectors
1351: Level: intermediate
1353: Note:
1354: `y` cannot be any of the `x` vectors.
1356: Developer Notes:
1357: This is a convenience routine, but implementations might be able to optimize it, for example, when `beta` is zero.
1359: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1360: @*/
1361: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
1362: {
1363: PetscFunctionBegin;
1365: VecCheckAssembled(y);
1367: PetscCall(VecSetErrorIfLocked(y, 1));
1368: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1371: if (y->ops->maxpby) {
1372: PetscInt zeros = 0;
1374: if (nv) {
1375: PetscAssertPointer(alpha, 3);
1376: PetscAssertPointer(x, 5);
1377: }
1379: for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
1383: PetscCheckSameTypeAndComm(y, 1, x[i], 5);
1384: VecCheckSameSize(y, 1, x[i], 5);
1385: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1386: VecCheckAssembled(x[i]);
1387: PetscCall(VecLockReadPush(x[i]));
1388: zeros += alpha[i] == (PetscScalar)0.0;
1389: }
1391: if (zeros < nv) { // has nonzero alpha
1392: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1393: PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1394: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1395: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1396: } else {
1397: PetscCall(VecScale(y, beta));
1398: }
1400: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1401: } else { // no maxpby
1402: if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1403: else PetscCall(VecScale(y, beta));
1404: PetscCall(VecMAXPY(y, nv, alpha, x));
1405: }
1406: PetscFunctionReturn(PETSC_SUCCESS);
1407: }
1409: /*@
1410: VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
1411: in the order they appear in the array. The concatenated vector resides on the same
1412: communicator and is the same type as the source vectors.
1414: Collective
1416: Input Parameters:
1417: + nx - number of vectors to be concatenated
1418: - X - array containing the vectors to be concatenated in the order of concatenation
1420: Output Parameters:
1421: + Y - concatenated vector
1422: - x_is - array of index sets corresponding to the concatenated components of `Y` (pass `NULL` if not needed)
1424: Level: advanced
1426: Notes:
1427: Concatenation is similar to the functionality of a `VECNEST` object; they both represent combination of
1428: different vector spaces. However, concatenated vectors do not store any information about their
1429: sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1430: manipulation of data in the concatenated vector that corresponds to the original components at creation.
1432: This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1433: has to operate on combined vector spaces and cannot utilize `VECNEST` objects due to incompatibility with
1434: bound projections.
1436: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1437: @*/
1438: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1439: {
1440: MPI_Comm comm;
1441: VecType vec_type;
1442: Vec Ytmp, Xtmp;
1443: IS *is_tmp;
1444: PetscInt i, shift = 0, Xnl, Xng, Xbegin;
1446: PetscFunctionBegin;
1450: PetscAssertPointer(Y, 3);
1452: if ((*X)->ops->concatenate) {
1453: /* use the dedicated concatenation function if available */
1454: PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1455: } else {
1456: /* loop over vectors and start creating IS */
1457: comm = PetscObjectComm((PetscObject)*X);
1458: PetscCall(VecGetType(*X, &vec_type));
1459: PetscCall(PetscMalloc1(nx, &is_tmp));
1460: for (i = 0; i < nx; i++) {
1461: PetscCall(VecGetSize(X[i], &Xng));
1462: PetscCall(VecGetLocalSize(X[i], &Xnl));
1463: PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1464: PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1465: shift += Xng;
1466: }
1467: /* create the concatenated vector */
1468: PetscCall(VecCreate(comm, &Ytmp));
1469: PetscCall(VecSetType(Ytmp, vec_type));
1470: PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1471: PetscCall(VecSetUp(Ytmp));
1472: /* copy data from X array to Y and return */
1473: for (i = 0; i < nx; i++) {
1474: PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1475: PetscCall(VecCopy(X[i], Xtmp));
1476: PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1477: }
1478: *Y = Ytmp;
1479: if (x_is) {
1480: *x_is = is_tmp;
1481: } else {
1482: for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1483: PetscCall(PetscFree(is_tmp));
1484: }
1485: }
1486: PetscFunctionReturn(PETSC_SUCCESS);
1487: }
1489: /* A helper function for VecGetSubVector to check if we can implement it with no-copy (i.e. the subvector shares
1490: memory with the original vector), and the block size of the subvector.
1492: Input Parameters:
1493: + X - the original vector
1494: - is - the index set of the subvector
1496: Output Parameters:
1497: + contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
1498: . start - start of contiguous block, as an offset from the start of the ownership range of the original vector
1499: - blocksize - the block size of the subvector
1501: */
1502: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1503: {
1504: PetscInt gstart, gend, lstart;
1505: PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1506: PetscInt n, N, ibs, vbs, bs = 1;
1508: PetscFunctionBegin;
1509: PetscCall(ISGetLocalSize(is, &n));
1510: PetscCall(ISGetSize(is, &N));
1511: PetscCall(ISGetBlockSize(is, &ibs));
1512: PetscCall(VecGetBlockSize(X, &vbs));
1513: PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1514: PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1515: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1516: if (ibs > 1) {
1517: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1518: bs = ibs;
1519: } else {
1520: if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1521: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1522: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1523: }
1525: *contig = red[0];
1526: *start = lstart;
1527: *blocksize = bs;
1528: PetscFunctionReturn(PETSC_SUCCESS);
1529: }
1531: /* A helper function for VecGetSubVector, to be used when we have to build a standalone subvector through VecScatter
1533: Input Parameters:
1534: + X - the original vector
1535: . is - the index set of the subvector
1536: - bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()
1538: Output Parameter:
1539: . Z - the subvector, which will compose the VecScatter context on output
1540: */
1541: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1542: {
1543: PetscInt n, N;
1544: VecScatter vscat;
1545: Vec Y;
1547: PetscFunctionBegin;
1548: PetscCall(ISGetLocalSize(is, &n));
1549: PetscCall(ISGetSize(is, &N));
1550: PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1551: PetscCall(VecSetSizes(Y, n, N));
1552: PetscCall(VecSetBlockSize(Y, bs));
1553: PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1554: PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1555: PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1556: PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1557: PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1558: PetscCall(VecScatterDestroy(&vscat));
1559: *Z = Y;
1560: PetscFunctionReturn(PETSC_SUCCESS);
1561: }
1563: /*@
1564: VecGetSubVector - Gets a vector representing part of another vector
1566: Collective
1568: Input Parameters:
1569: + X - vector from which to extract a subvector
1570: - is - index set representing portion of `X` to extract
1572: Output Parameter:
1573: . Y - subvector corresponding to `is`
1575: Level: advanced
1577: Notes:
1578: The subvector `Y` should be returned with `VecRestoreSubVector()`.
1579: `X` and `is` must be defined on the same communicator
1581: Changes to the subvector will be reflected in the `X` vector on the call to `VecRestoreSubVector()`.
1583: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1584: modifying the subvector. Other non-overlapping subvectors can still be obtained from `X` using this function.
1586: 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`.
1588: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1589: @*/
1590: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1591: {
1592: Vec Z;
1594: PetscFunctionBegin;
1597: PetscCheckSameComm(X, 1, is, 2);
1598: PetscAssertPointer(Y, 3);
1599: if (X->ops->getsubvector) {
1600: PetscUseTypeMethod(X, getsubvector, is, &Z);
1601: } else { /* Default implementation currently does no caching */
1602: PetscBool contig;
1603: PetscInt n, N, start, bs;
1605: PetscCall(ISGetLocalSize(is, &n));
1606: PetscCall(ISGetSize(is, &N));
1607: PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1608: if (contig) { /* We can do a no-copy implementation */
1609: const PetscScalar *x;
1610: PetscInt state = 0;
1611: PetscBool isstd, iscuda, iship;
1613: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1614: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1615: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1616: if (iscuda) {
1617: #if defined(PETSC_HAVE_CUDA)
1618: const PetscScalar *x_d;
1619: PetscMPIInt size;
1620: PetscOffloadMask flg;
1622: PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1623: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1624: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1625: if (x) x += start;
1626: if (x_d) x_d += start;
1627: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1628: if (size == 1) {
1629: PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1630: } else {
1631: PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1632: }
1633: Z->offloadmask = flg;
1634: #endif
1635: } else if (iship) {
1636: #if defined(PETSC_HAVE_HIP)
1637: const PetscScalar *x_d;
1638: PetscMPIInt size;
1639: PetscOffloadMask flg;
1641: PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1642: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1643: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1644: if (x) x += start;
1645: if (x_d) x_d += start;
1646: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1647: if (size == 1) {
1648: PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1649: } else {
1650: PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1651: }
1652: Z->offloadmask = flg;
1653: #endif
1654: } else if (isstd) {
1655: PetscMPIInt size;
1657: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1658: PetscCall(VecGetArrayRead(X, &x));
1659: if (x) x += start;
1660: if (size == 1) {
1661: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1662: } else {
1663: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1664: }
1665: PetscCall(VecRestoreArrayRead(X, &x));
1666: } else { /* default implementation: use place array */
1667: PetscCall(VecGetArrayRead(X, &x));
1668: PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1669: PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1670: PetscCall(VecSetSizes(Z, n, N));
1671: PetscCall(VecSetBlockSize(Z, bs));
1672: PetscCall(VecPlaceArray(Z, PetscSafePointerPlusOffset(x, start)));
1673: PetscCall(VecRestoreArrayRead(X, &x));
1674: }
1676: /* this is relevant only in debug mode */
1677: PetscCall(VecLockGet(X, &state));
1678: if (state) PetscCall(VecLockReadPush(Z));
1679: Z->ops->placearray = NULL;
1680: Z->ops->replacearray = NULL;
1681: } else { /* Have to create a scatter and do a copy */
1682: PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1683: }
1684: }
1685: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1686: if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1687: PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1688: *Y = Z;
1689: PetscFunctionReturn(PETSC_SUCCESS);
1690: }
1692: /*@
1693: VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`
1695: Collective
1697: Input Parameters:
1698: + X - vector from which subvector was obtained
1699: . is - index set representing the subset of `X`
1700: - Y - subvector being restored
1702: Level: advanced
1704: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1705: @*/
1706: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1707: {
1708: PETSC_UNUSED PetscObjectState dummystate = 0;
1709: PetscBool unchanged;
1711: PetscFunctionBegin;
1714: PetscCheckSameComm(X, 1, is, 2);
1715: PetscAssertPointer(Y, 3);
1718: if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1719: else {
1720: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1721: if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1722: VecScatter scatter;
1723: PetscInt state;
1725: PetscCall(VecLockGet(X, &state));
1726: PetscCheck(state == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec X is locked for read-only or read/write access");
1728: PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1729: if (scatter) {
1730: PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1731: PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1732: } else {
1733: PetscBool iscuda, iship;
1734: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1735: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1737: if (iscuda) {
1738: #if defined(PETSC_HAVE_CUDA)
1739: PetscOffloadMask ymask = (*Y)->offloadmask;
1741: /* The offloadmask of X dictates where to move memory
1742: If X GPU data is valid, then move Y data on GPU if needed
1743: Otherwise, move back to the CPU */
1744: switch (X->offloadmask) {
1745: case PETSC_OFFLOAD_BOTH:
1746: if (ymask == PETSC_OFFLOAD_CPU) {
1747: PetscCall(VecCUDAResetArray(*Y));
1748: } else if (ymask == PETSC_OFFLOAD_GPU) {
1749: X->offloadmask = PETSC_OFFLOAD_GPU;
1750: }
1751: break;
1752: case PETSC_OFFLOAD_GPU:
1753: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1754: break;
1755: case PETSC_OFFLOAD_CPU:
1756: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1757: break;
1758: case PETSC_OFFLOAD_UNALLOCATED:
1759: case PETSC_OFFLOAD_KOKKOS:
1760: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1761: }
1762: #endif
1763: } else if (iship) {
1764: #if defined(PETSC_HAVE_HIP)
1765: PetscOffloadMask ymask = (*Y)->offloadmask;
1767: /* The offloadmask of X dictates where to move memory
1768: If X GPU data is valid, then move Y data on GPU if needed
1769: Otherwise, move back to the CPU */
1770: switch (X->offloadmask) {
1771: case PETSC_OFFLOAD_BOTH:
1772: if (ymask == PETSC_OFFLOAD_CPU) {
1773: PetscCall(VecHIPResetArray(*Y));
1774: } else if (ymask == PETSC_OFFLOAD_GPU) {
1775: X->offloadmask = PETSC_OFFLOAD_GPU;
1776: }
1777: break;
1778: case PETSC_OFFLOAD_GPU:
1779: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1780: break;
1781: case PETSC_OFFLOAD_CPU:
1782: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1783: break;
1784: case PETSC_OFFLOAD_UNALLOCATED:
1785: case PETSC_OFFLOAD_KOKKOS:
1786: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1787: }
1788: #endif
1789: } else {
1790: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1791: PetscCall(VecResetArray(*Y));
1792: }
1793: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1794: }
1795: }
1796: }
1797: PetscCall(VecDestroy(Y));
1798: PetscFunctionReturn(PETSC_SUCCESS);
1799: }
1801: /*@
1802: VecCreateLocalVector - Creates a vector object suitable for use with `VecGetLocalVector()` and friends. You must call `VecDestroy()` when the
1803: vector is no longer needed.
1805: Not Collective.
1807: Input Parameter:
1808: . v - The vector for which the local vector is desired.
1810: Output Parameter:
1811: . w - Upon exit this contains the local vector.
1813: Level: beginner
1815: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1816: @*/
1817: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1818: {
1819: VecType roottype;
1820: PetscInt n;
1822: PetscFunctionBegin;
1824: PetscAssertPointer(w, 2);
1825: if (v->ops->createlocalvector) {
1826: PetscUseTypeMethod(v, createlocalvector, w);
1827: PetscFunctionReturn(PETSC_SUCCESS);
1828: }
1829: PetscCall(VecGetRootType_Private(v, &roottype));
1830: PetscCall(VecCreate(PETSC_COMM_SELF, w));
1831: PetscCall(VecGetLocalSize(v, &n));
1832: PetscCall(VecSetSizes(*w, n, n));
1833: PetscCall(VecGetBlockSize(v, &n));
1834: PetscCall(VecSetBlockSize(*w, n));
1835: PetscCall(VecSetType(*w, roottype));
1836: PetscFunctionReturn(PETSC_SUCCESS);
1837: }
1839: /*@
1840: VecGetLocalVectorRead - Maps the local portion of a vector into a
1841: vector.
1843: Not Collective.
1845: Input Parameter:
1846: . v - The vector for which the local vector is desired.
1848: Output Parameter:
1849: . w - Upon exit this contains the local vector.
1851: Level: beginner
1853: Notes:
1854: You must call `VecRestoreLocalVectorRead()` when the local
1855: vector is no longer needed.
1857: This function is similar to `VecGetArrayRead()` which maps the local
1858: portion into a raw pointer. `VecGetLocalVectorRead()` is usually
1859: almost as efficient as `VecGetArrayRead()` but in certain circumstances
1860: `VecGetLocalVectorRead()` can be much more efficient than
1861: `VecGetArrayRead()`. This is because the construction of a contiguous
1862: array representing the vector data required by `VecGetArrayRead()` can
1863: be an expensive operation for certain vector types. For example, for
1864: GPU vectors `VecGetArrayRead()` requires that the data between device
1865: and host is synchronized.
1867: Unlike `VecGetLocalVector()`, this routine is not collective and
1868: preserves cached information.
1870: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1871: @*/
1872: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1873: {
1874: PetscFunctionBegin;
1877: VecCheckSameLocalSize(v, 1, w, 2);
1878: if (v->ops->getlocalvectorread) {
1879: PetscUseTypeMethod(v, getlocalvectorread, w);
1880: } else {
1881: PetscScalar *a;
1883: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1884: PetscCall(VecPlaceArray(w, a));
1885: }
1886: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1887: PetscCall(VecLockReadPush(v));
1888: PetscCall(VecLockReadPush(w));
1889: PetscFunctionReturn(PETSC_SUCCESS);
1890: }
1892: /*@
1893: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1894: previously mapped into a vector using `VecGetLocalVectorRead()`.
1896: Not Collective.
1898: Input Parameters:
1899: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
1900: - w - The vector into which the local portion of `v` was mapped.
1902: Level: beginner
1904: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1905: @*/
1906: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1907: {
1908: PetscFunctionBegin;
1911: if (v->ops->restorelocalvectorread) {
1912: PetscUseTypeMethod(v, restorelocalvectorread, w);
1913: } else {
1914: const PetscScalar *a;
1916: PetscCall(VecGetArrayRead(w, &a));
1917: PetscCall(VecRestoreArrayRead(v, &a));
1918: PetscCall(VecResetArray(w));
1919: }
1920: PetscCall(VecLockReadPop(v));
1921: PetscCall(VecLockReadPop(w));
1922: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1923: PetscFunctionReturn(PETSC_SUCCESS);
1924: }
1926: /*@
1927: VecGetLocalVector - Maps the local portion of a vector into a
1928: vector.
1930: Collective
1932: Input Parameter:
1933: . v - The vector for which the local vector is desired.
1935: Output Parameter:
1936: . w - Upon exit this contains the local vector.
1938: Level: beginner
1940: Notes:
1941: You must call `VecRestoreLocalVector()` when the local
1942: vector is no longer needed.
1944: This function is similar to `VecGetArray()` which maps the local
1945: portion into a raw pointer. `VecGetLocalVector()` is usually about as
1946: efficient as `VecGetArray()` but in certain circumstances
1947: `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1948: This is because the construction of a contiguous array representing
1949: the vector data required by `VecGetArray()` can be an expensive
1950: operation for certain vector types. For example, for GPU vectors
1951: `VecGetArray()` requires that the data between device and host is
1952: synchronized.
1954: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1955: @*/
1956: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1957: {
1958: PetscFunctionBegin;
1961: VecCheckSameLocalSize(v, 1, w, 2);
1962: if (v->ops->getlocalvector) {
1963: PetscUseTypeMethod(v, getlocalvector, w);
1964: } else {
1965: PetscScalar *a;
1967: PetscCall(VecGetArray(v, &a));
1968: PetscCall(VecPlaceArray(w, a));
1969: }
1970: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1971: PetscFunctionReturn(PETSC_SUCCESS);
1972: }
1974: /*@
1975: VecRestoreLocalVector - Unmaps the local portion of a vector
1976: previously mapped into a vector using `VecGetLocalVector()`.
1978: Logically Collective.
1980: Input Parameters:
1981: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
1982: - w - The vector into which the local portion of `v` was mapped.
1984: Level: beginner
1986: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1987: @*/
1988: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1989: {
1990: PetscFunctionBegin;
1993: if (v->ops->restorelocalvector) {
1994: PetscUseTypeMethod(v, restorelocalvector, w);
1995: } else {
1996: PetscScalar *a;
1997: PetscCall(VecGetArray(w, &a));
1998: PetscCall(VecRestoreArray(v, &a));
1999: PetscCall(VecResetArray(w));
2000: }
2001: PetscCall(PetscObjectStateIncrease((PetscObject)w));
2002: PetscCall(PetscObjectStateIncrease((PetscObject)v));
2003: PetscFunctionReturn(PETSC_SUCCESS);
2004: }
2006: /*@C
2007: VecGetArray - Returns a pointer to a contiguous array that contains this
2008: MPI processes's portion of the vector data
2010: Logically Collective
2012: Input Parameter:
2013: . x - the vector
2015: Output Parameter:
2016: . a - location to put pointer to the array
2018: Level: beginner
2020: Notes:
2021: For the standard PETSc vectors, `VecGetArray()` returns a pointer to the local data array and
2022: does not use any copies. If the underlying vector data is not stored in a contiguous array
2023: this routine will copy the data to a contiguous array and return a pointer to that. You MUST
2024: call `VecRestoreArray()` when you no longer need access to the array.
2026: For vectors that may also have the array data in GPU memory, for example, `VECCUDA`, this call ensures the CPU array has the
2027: most recent array values by copying the data from the GPU memory if needed.
2029: Fortran Note:
2030: .vb
2031: PetscScalar, pointer :: a(:)
2032: .ve
2034: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecPlaceArray()`, `VecGetArray2d()`,
2035: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`, `VecGetArrayAndMemType()`
2036: @*/
2037: PetscErrorCode VecGetArray(Vec x, PetscScalar *a[])
2038: {
2039: PetscFunctionBegin;
2041: PetscCall(VecSetErrorIfLocked(x, 1));
2042: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
2043: PetscUseTypeMethod(x, getarray, a);
2044: } else if (x->petscnative) { /* VECSTANDARD */
2045: *a = *((PetscScalar **)x->data);
2046: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
2047: PetscFunctionReturn(PETSC_SUCCESS);
2048: }
2050: /*@C
2051: VecRestoreArray - Restores a vector after `VecGetArray()` has been called and the array is no longer needed
2053: Logically Collective
2055: Input Parameters:
2056: + x - the vector
2057: - a - location of pointer to array obtained from `VecGetArray()`
2059: Level: beginner
2061: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2062: `VecGetArrayPair()`, `VecRestoreArrayPair()`
2063: @*/
2064: PetscErrorCode VecRestoreArray(Vec x, PetscScalar *a[])
2065: {
2066: PetscFunctionBegin;
2068: if (a) PetscAssertPointer(a, 2);
2069: if (x->ops->restorearray) {
2070: PetscUseTypeMethod(x, restorearray, a);
2071: } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
2072: if (a) *a = NULL;
2073: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2074: PetscFunctionReturn(PETSC_SUCCESS);
2075: }
2076: /*@C
2077: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
2079: Not Collective
2081: Input Parameter:
2082: . x - the vector
2084: Output Parameter:
2085: . a - the array
2087: Level: beginner
2089: Notes:
2090: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2092: Unlike `VecGetArray()`, preserves cached information like vector norms.
2094: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
2095: implementations may require a copy, but such implementations should cache the contiguous representation so that
2096: only one copy is performed when this routine is called multiple times in sequence.
2098: For vectors that may also have the array data in GPU memory, for example, `VECCUDA`, this call ensures the CPU array has the
2099: most recent array values by copying the data from the GPU memory if needed.
2101: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2102: `VecGetArrayAndMemType()`
2103: @*/
2104: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar *a[])
2105: {
2106: PetscFunctionBegin;
2108: PetscAssertPointer(a, 2);
2109: if (x->ops->getarrayread) {
2110: PetscUseTypeMethod(x, getarrayread, a);
2111: } else if (x->ops->getarray) {
2112: PetscObjectState state;
2114: /* VECNEST, VECCUDA, VECKOKKOS etc */
2115: // x->ops->getarray may bump the object state, but since we know this is a read-only get
2116: // we can just undo that
2117: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2118: PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
2119: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2120: } else if (x->petscnative) {
2121: /* VECSTANDARD */
2122: *a = *((PetscScalar **)x->data);
2123: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
2124: PetscFunctionReturn(PETSC_SUCCESS);
2125: }
2127: /*@C
2128: VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`
2130: Not Collective
2132: Input Parameters:
2133: + x - the vector
2134: - a - the array
2136: Level: beginner
2138: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2139: @*/
2140: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar *a[])
2141: {
2142: PetscFunctionBegin;
2144: if (a) PetscAssertPointer(a, 2);
2145: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
2146: /* nothing */
2147: } else if (x->ops->restorearrayread) { /* VECNEST */
2148: PetscUseTypeMethod(x, restorearrayread, a);
2149: } else { /* No one? */
2150: PetscObjectState state;
2152: // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2153: // we can just undo that
2154: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2155: PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2156: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2157: }
2158: if (a) *a = NULL;
2159: PetscFunctionReturn(PETSC_SUCCESS);
2160: }
2162: /*@C
2163: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
2164: MPI processes's portion of the vector data.
2166: Logically Collective
2168: Input Parameter:
2169: . x - the vector
2171: Output Parameter:
2172: . a - location to put pointer to the array
2174: Level: intermediate
2176: Note:
2177: The values in this array are NOT valid, the caller of this routine is responsible for putting
2178: values into the array; any values it does not set will be invalid.
2180: The array must be returned using a matching call to `VecRestoreArrayWrite()`.
2182: For vectors associated with GPUs, the host and device vectors are not synchronized before
2183: giving access. If you need correct values in the array use `VecGetArray()`
2185: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecPlaceArray()`, `VecGetArray2d()`,
2186: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`, `VecGetArrayAndMemType()`
2187: @*/
2188: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar *a[])
2189: {
2190: PetscFunctionBegin;
2192: PetscAssertPointer(a, 2);
2193: PetscCall(VecSetErrorIfLocked(x, 1));
2194: if (x->ops->getarraywrite) {
2195: PetscUseTypeMethod(x, getarraywrite, a);
2196: } else {
2197: PetscCall(VecGetArray(x, a));
2198: }
2199: PetscFunctionReturn(PETSC_SUCCESS);
2200: }
2202: /*@C
2203: VecRestoreArrayWrite - Restores a vector after `VecGetArrayWrite()` has been called.
2205: Logically Collective
2207: Input Parameters:
2208: + x - the vector
2209: - a - location of pointer to array obtained from `VecGetArray()`
2211: Level: beginner
2213: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2214: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2215: @*/
2216: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar *a[])
2217: {
2218: PetscFunctionBegin;
2220: if (a) PetscAssertPointer(a, 2);
2221: if (x->ops->restorearraywrite) {
2222: PetscUseTypeMethod(x, restorearraywrite, a);
2223: } else if (x->ops->restorearray) {
2224: PetscUseTypeMethod(x, restorearray, a);
2225: }
2226: if (a) *a = NULL;
2227: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2228: PetscFunctionReturn(PETSC_SUCCESS);
2229: }
2231: /*@C
2232: VecGetArrays - Returns a pointer to the arrays in a set of vectors
2233: that were created by a call to `VecDuplicateVecs()`.
2235: Logically Collective; No Fortran Support
2237: Input Parameters:
2238: + x - the vectors
2239: - n - the number of vectors
2241: Output Parameter:
2242: . a - location to put pointer to the array
2244: Level: intermediate
2246: Note:
2247: You MUST call `VecRestoreArrays()` when you no longer need access to the arrays.
2249: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2250: @*/
2251: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2252: {
2253: PetscInt i;
2254: PetscScalar **q;
2256: PetscFunctionBegin;
2257: PetscAssertPointer(x, 1);
2259: PetscAssertPointer(a, 3);
2260: PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2261: PetscCall(PetscMalloc1(n, &q));
2262: for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2263: *a = q;
2264: PetscFunctionReturn(PETSC_SUCCESS);
2265: }
2267: /*@C
2268: VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2269: has been called.
2271: Logically Collective; No Fortran Support
2273: Input Parameters:
2274: + x - the vector
2275: . n - the number of vectors
2276: - a - location of pointer to arrays obtained from `VecGetArrays()`
2278: Notes:
2279: For regular PETSc vectors this routine does not involve any copies. For
2280: any special vectors that do not store local vector data in a contiguous
2281: array, this routine will copy the data back into the underlying
2282: vector data structure from the arrays obtained with `VecGetArrays()`.
2284: Level: intermediate
2286: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2287: @*/
2288: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2289: {
2290: PetscInt i;
2291: PetscScalar **q = *a;
2293: PetscFunctionBegin;
2294: PetscAssertPointer(x, 1);
2296: PetscAssertPointer(a, 3);
2298: for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2299: PetscCall(PetscFree(q));
2300: PetscFunctionReturn(PETSC_SUCCESS);
2301: }
2303: /*@C
2304: VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g.,
2305: `VECCUDA`), the returned pointer will be a device pointer to the device memory that contains
2306: this MPI processes's portion of the vector data.
2308: Logically Collective; No Fortran Support
2310: Input Parameter:
2311: . x - the vector
2313: Output Parameters:
2314: + a - location to put pointer to the array
2315: - mtype - memory type of the array
2317: Level: beginner
2319: Note:
2320: Device data is guaranteed to have the latest value. Otherwise, when this is a host vector
2321: (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host
2322: pointer.
2324: For `VECKOKKOS`, if Kokkos is configured without device (e.g., use serial or openmp), per
2325: this function, the vector works like `VECSEQ`/`VECMPI`; otherwise, it works like `VECCUDA` or
2326: `VECHIP` etc.
2328: Use `VecRestoreArrayAndMemType()` when the array access is no longer needed.
2330: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`,
2331: `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2332: @*/
2333: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar *a[], PetscMemType *mtype)
2334: {
2335: PetscFunctionBegin;
2338: if (a) PetscAssertPointer(a, 2);
2339: if (mtype) PetscAssertPointer(mtype, 3);
2340: PetscCall(VecSetErrorIfLocked(x, 1));
2341: if (x->ops->getarrayandmemtype) {
2342: /* VECCUDA, VECKOKKOS etc */
2343: PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2344: } else {
2345: /* VECSTANDARD, VECNEST, VECVIENNACL */
2346: PetscCall(VecGetArray(x, a));
2347: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2348: }
2349: PetscFunctionReturn(PETSC_SUCCESS);
2350: }
2352: /*@C
2353: VecRestoreArrayAndMemType - Restores a vector after `VecGetArrayAndMemType()` has been called.
2355: Logically Collective; No Fortran Support
2357: Input Parameters:
2358: + x - the vector
2359: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`
2361: Level: beginner
2363: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`,
2364: `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2365: @*/
2366: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar *a[])
2367: {
2368: PetscFunctionBegin;
2371: if (a) PetscAssertPointer(a, 2);
2372: if (x->ops->restorearrayandmemtype) {
2373: /* VECCUDA, VECKOKKOS etc */
2374: PetscUseTypeMethod(x, restorearrayandmemtype, a);
2375: } else {
2376: /* VECNEST, VECVIENNACL */
2377: PetscCall(VecRestoreArray(x, a));
2378: } /* VECSTANDARD does nothing */
2379: if (a) *a = NULL;
2380: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2381: PetscFunctionReturn(PETSC_SUCCESS);
2382: }
2384: /*@C
2385: VecGetArrayReadAndMemType - Like `VecGetArrayRead()`, but if the input vector is a device vector, it will return a read-only device pointer.
2386: The returned pointer is guaranteed to point to up-to-date data. For host vectors, it functions as `VecGetArrayRead()`.
2388: Not Collective; No Fortran Support
2390: Input Parameter:
2391: . x - the vector
2393: Output Parameters:
2394: + a - the array
2395: - mtype - memory type of the array
2397: Level: beginner
2399: Notes:
2400: The array must be returned using a matching call to `VecRestoreArrayReadAndMemType()`.
2402: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2403: @*/
2404: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar *a[], PetscMemType *mtype)
2405: {
2406: PetscFunctionBegin;
2409: PetscAssertPointer(a, 2);
2410: if (mtype) PetscAssertPointer(mtype, 3);
2411: if (x->ops->getarrayreadandmemtype) {
2412: /* VECCUDA/VECHIP though they are also petscnative */
2413: PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2414: } else if (x->ops->getarrayandmemtype) {
2415: /* VECKOKKOS */
2416: PetscObjectState state;
2418: // see VecGetArrayRead() for why
2419: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2420: PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2421: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2422: } else {
2423: PetscCall(VecGetArrayRead(x, a));
2424: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2425: }
2426: PetscFunctionReturn(PETSC_SUCCESS);
2427: }
2429: /*@C
2430: VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`
2432: Not Collective; No Fortran Support
2434: Input Parameters:
2435: + x - the vector
2436: - a - the array
2438: Level: beginner
2440: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2441: @*/
2442: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar *a[])
2443: {
2444: PetscFunctionBegin;
2447: if (a) PetscAssertPointer(a, 2);
2448: if (x->ops->restorearrayreadandmemtype) {
2449: /* VECCUDA/VECHIP */
2450: PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2451: } else if (!x->petscnative) {
2452: /* VECNEST */
2453: PetscCall(VecRestoreArrayRead(x, a));
2454: }
2455: if (a) *a = NULL;
2456: PetscFunctionReturn(PETSC_SUCCESS);
2457: }
2459: /*@C
2460: VecGetArrayWriteAndMemType - Like `VecGetArrayWrite()`, but if this is a device vector it will always return
2461: a device pointer to the device memory that contains this processor's portion of the vector data.
2463: Logically Collective; No Fortran Support
2465: Input Parameter:
2466: . x - the vector
2468: Output Parameters:
2469: + a - the array
2470: - mtype - memory type of the array
2472: Level: beginner
2474: Note:
2475: The array must be returned using a matching call to `VecRestoreArrayWriteAndMemType()`, where it will label the device memory as most recent.
2477: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2478: @*/
2479: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar *a[], PetscMemType *mtype)
2480: {
2481: PetscFunctionBegin;
2484: PetscCall(VecSetErrorIfLocked(x, 1));
2485: PetscAssertPointer(a, 2);
2486: if (mtype) PetscAssertPointer(mtype, 3);
2487: if (x->ops->getarraywriteandmemtype) {
2488: /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2489: PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2490: } else if (x->ops->getarrayandmemtype) {
2491: PetscCall(VecGetArrayAndMemType(x, a, mtype));
2492: } else {
2493: /* VECNEST, VECVIENNACL */
2494: PetscCall(VecGetArrayWrite(x, a));
2495: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2496: }
2497: PetscFunctionReturn(PETSC_SUCCESS);
2498: }
2500: /*@C
2501: VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`
2503: Logically Collective; No Fortran Support
2505: Input Parameters:
2506: + x - the vector
2507: - a - the array
2509: Level: beginner
2511: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2512: @*/
2513: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar *a[])
2514: {
2515: PetscFunctionBegin;
2518: PetscCall(VecSetErrorIfLocked(x, 1));
2519: if (a) PetscAssertPointer(a, 2);
2520: if (x->ops->restorearraywriteandmemtype) {
2521: /* VECCUDA/VECHIP */
2522: PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2523: PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2524: } else if (x->ops->restorearrayandmemtype) {
2525: PetscCall(VecRestoreArrayAndMemType(x, a));
2526: } else {
2527: PetscCall(VecRestoreArray(x, a));
2528: }
2529: if (a) *a = NULL;
2530: PetscFunctionReturn(PETSC_SUCCESS);
2531: }
2533: /*@
2534: VecPlaceArray - Allows one to replace the array in a vector with an
2535: array provided by the user. This is useful to avoid copying an array
2536: into a vector.
2538: Logically Collective; No Fortran Support
2540: Input Parameters:
2541: + vec - the vector
2542: - array - the array
2544: Level: developer
2546: Notes:
2547: Adding `const` to `array` was an oversight, as subsequent operations on `vec` would
2548: likely modify the data in `array`. However, we have kept it to avoid breaking APIs.
2550: Use `VecReplaceArray()` instead to permanently replace the array
2552: You can return to the original array with a call to `VecResetArray()`. `vec` does not take
2553: ownership of `array` in any way.
2555: The user must free `array` themselves but be careful not to
2556: do so before the vector has either been destroyed, had its original array restored with
2557: `VecResetArray()` or permanently replaced with `VecReplaceArray()`.
2559: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2560: @*/
2561: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2562: {
2563: PetscFunctionBegin;
2566: if (array) PetscAssertPointer(array, 2);
2567: PetscUseTypeMethod(vec, placearray, array);
2568: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2569: PetscFunctionReturn(PETSC_SUCCESS);
2570: }
2572: /*@C
2573: VecReplaceArray - Allows one to replace the array in a vector with an
2574: array provided by the user. This is useful to avoid copying an array
2575: into a vector.
2577: Logically Collective; No Fortran Support
2579: Input Parameters:
2580: + vec - the vector
2581: - array - the array
2583: Level: developer
2585: Notes:
2586: Adding `const` to `array` was an oversight, as subsequent operations on `vec` would
2587: likely modify the data in `array`. However, we have kept it to avoid breaking APIs.
2589: This permanently replaces the array and frees the memory associated
2590: with the old array. Use `VecPlaceArray()` to temporarily replace the array.
2592: The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
2593: freed by the user. It will be freed when the vector is destroyed.
2595: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2596: @*/
2597: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2598: {
2599: PetscFunctionBegin;
2602: PetscUseTypeMethod(vec, replacearray, array);
2603: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2604: PetscFunctionReturn(PETSC_SUCCESS);
2605: }
2607: /*@C
2608: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2609: processor's portion of the vector data. You MUST call `VecRestoreArray2d()`
2610: when you no longer need access to the array.
2612: Logically Collective
2614: Input Parameters:
2615: + x - the vector
2616: . m - first dimension of two dimensional array
2617: . n - second dimension of two dimensional array
2618: . mstart - first index you will use in first coordinate direction (often 0)
2619: - nstart - first index in the second coordinate direction (often 0)
2621: Output Parameter:
2622: . a - location to put pointer to the array
2624: Level: developer
2626: Notes:
2627: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2628: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2629: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2630: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2632: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2634: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2635: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2636: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2637: @*/
2638: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2639: {
2640: PetscInt i, N;
2641: PetscScalar *aa;
2643: PetscFunctionBegin;
2645: PetscAssertPointer(a, 6);
2647: PetscCall(VecGetLocalSize(x, &N));
2648: 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);
2649: PetscCall(VecGetArray(x, &aa));
2651: PetscCall(PetscMalloc1(m, a));
2652: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2653: *a -= mstart;
2654: PetscFunctionReturn(PETSC_SUCCESS);
2655: }
2657: /*@C
2658: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2659: processor's portion of the vector data. You MUST call `VecRestoreArray2dWrite()`
2660: when you no longer need access to the array.
2662: Logically Collective
2664: Input Parameters:
2665: + x - the vector
2666: . m - first dimension of two dimensional array
2667: . n - second dimension of two dimensional array
2668: . mstart - first index you will use in first coordinate direction (often 0)
2669: - nstart - first index in the second coordinate direction (often 0)
2671: Output Parameter:
2672: . a - location to put pointer to the array
2674: Level: developer
2676: Notes:
2677: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2678: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2679: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2680: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2682: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2684: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2685: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2686: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2687: @*/
2688: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2689: {
2690: PetscInt i, N;
2691: PetscScalar *aa;
2693: PetscFunctionBegin;
2695: PetscAssertPointer(a, 6);
2697: PetscCall(VecGetLocalSize(x, &N));
2698: 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);
2699: PetscCall(VecGetArrayWrite(x, &aa));
2701: PetscCall(PetscMalloc1(m, a));
2702: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2703: *a -= mstart;
2704: PetscFunctionReturn(PETSC_SUCCESS);
2705: }
2707: /*@C
2708: VecRestoreArray2d - Restores a vector after `VecGetArray2d()` has been called.
2710: Logically Collective
2712: Input Parameters:
2713: + x - the vector
2714: . m - first dimension of two dimensional array
2715: . n - second dimension of the two dimensional array
2716: . mstart - first index you will use in first coordinate direction (often 0)
2717: . nstart - first index in the second coordinate direction (often 0)
2718: - a - location of pointer to array obtained from `VecGetArray2d()`
2720: Level: developer
2722: Notes:
2723: For regular PETSc vectors this routine does not involve any copies. For
2724: any special vectors that do not store local vector data in a contiguous
2725: array, this routine will copy the data back into the underlying
2726: vector data structure from the array obtained with `VecGetArray()`.
2728: This routine actually zeros out the `a` pointer.
2730: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2731: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2732: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2733: @*/
2734: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2735: {
2736: void *dummy;
2738: PetscFunctionBegin;
2740: PetscAssertPointer(a, 6);
2742: dummy = (void *)(*a + mstart);
2743: PetscCall(PetscFree(dummy));
2744: PetscCall(VecRestoreArray(x, NULL));
2745: *a = NULL;
2746: PetscFunctionReturn(PETSC_SUCCESS);
2747: }
2749: /*@C
2750: VecRestoreArray2dWrite - Restores a vector after `VecGetArray2dWrite()` has been called.
2752: Logically Collective
2754: Input Parameters:
2755: + x - the vector
2756: . m - first dimension of two dimensional array
2757: . n - second dimension of the two dimensional array
2758: . mstart - first index you will use in first coordinate direction (often 0)
2759: . nstart - first index in the second coordinate direction (often 0)
2760: - a - location of pointer to array obtained from `VecGetArray2d()`
2762: Level: developer
2764: Notes:
2765: For regular PETSc vectors this routine does not involve any copies. For
2766: any special vectors that do not store local vector data in a contiguous
2767: array, this routine will copy the data back into the underlying
2768: vector data structure from the array obtained with `VecGetArray()`.
2770: This routine actually zeros out the `a` pointer.
2772: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2773: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2774: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2775: @*/
2776: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2777: {
2778: void *dummy;
2780: PetscFunctionBegin;
2782: PetscAssertPointer(a, 6);
2784: dummy = (void *)(*a + mstart);
2785: PetscCall(PetscFree(dummy));
2786: PetscCall(VecRestoreArrayWrite(x, NULL));
2787: PetscFunctionReturn(PETSC_SUCCESS);
2788: }
2790: /*@C
2791: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2792: processor's portion of the vector data. You MUST call `VecRestoreArray1d()`
2793: when you no longer need access to the array.
2795: Logically Collective
2797: Input Parameters:
2798: + x - the vector
2799: . m - first dimension of two dimensional array
2800: - mstart - first index you will use in first coordinate direction (often 0)
2802: Output Parameter:
2803: . a - location to put pointer to the array
2805: Level: developer
2807: Notes:
2808: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2809: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2810: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2812: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2814: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2815: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2816: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2817: @*/
2818: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2819: {
2820: PetscInt N;
2822: PetscFunctionBegin;
2824: PetscAssertPointer(a, 4);
2826: PetscCall(VecGetLocalSize(x, &N));
2827: 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);
2828: PetscCall(VecGetArray(x, a));
2829: *a -= mstart;
2830: PetscFunctionReturn(PETSC_SUCCESS);
2831: }
2833: /*@C
2834: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
2835: processor's portion of the vector data. You MUST call `VecRestoreArray1dWrite()`
2836: when you no longer need access to the array.
2838: Logically Collective
2840: Input Parameters:
2841: + x - the vector
2842: . m - first dimension of two dimensional array
2843: - mstart - first index you will use in first coordinate direction (often 0)
2845: Output Parameter:
2846: . a - location to put pointer to the array
2848: Level: developer
2850: Notes:
2851: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2852: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2853: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2855: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2857: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2858: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2859: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2860: @*/
2861: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2862: {
2863: PetscInt N;
2865: PetscFunctionBegin;
2867: PetscAssertPointer(a, 4);
2869: PetscCall(VecGetLocalSize(x, &N));
2870: 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);
2871: PetscCall(VecGetArrayWrite(x, a));
2872: *a -= mstart;
2873: PetscFunctionReturn(PETSC_SUCCESS);
2874: }
2876: /*@C
2877: VecRestoreArray1d - Restores a vector after `VecGetArray1d()` has been called.
2879: Logically Collective
2881: Input Parameters:
2882: + x - the vector
2883: . m - first dimension of two dimensional array
2884: . mstart - first index you will use in first coordinate direction (often 0)
2885: - a - location of pointer to array obtained from `VecGetArray1d()`
2887: Level: developer
2889: Notes:
2890: For regular PETSc vectors this routine does not involve any copies. For
2891: any special vectors that do not store local vector data in a contiguous
2892: array, this routine will copy the data back into the underlying
2893: vector data structure from the array obtained with `VecGetArray1d()`.
2895: This routine actually zeros out the `a` pointer.
2897: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2898: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2899: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2900: @*/
2901: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2902: {
2903: PetscFunctionBegin;
2906: PetscCall(VecRestoreArray(x, NULL));
2907: *a = NULL;
2908: PetscFunctionReturn(PETSC_SUCCESS);
2909: }
2911: /*@C
2912: VecRestoreArray1dWrite - Restores a vector after `VecGetArray1dWrite()` has been called.
2914: Logically Collective
2916: Input Parameters:
2917: + x - the vector
2918: . m - first dimension of two dimensional array
2919: . mstart - first index you will use in first coordinate direction (often 0)
2920: - a - location of pointer to array obtained from `VecGetArray1d()`
2922: Level: developer
2924: Notes:
2925: For regular PETSc vectors this routine does not involve any copies. For
2926: any special vectors that do not store local vector data in a contiguous
2927: array, this routine will copy the data back into the underlying
2928: vector data structure from the array obtained with `VecGetArray1d()`.
2930: This routine actually zeros out the `a` pointer.
2932: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2933: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2934: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2935: @*/
2936: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2937: {
2938: PetscFunctionBegin;
2941: PetscCall(VecRestoreArrayWrite(x, NULL));
2942: *a = NULL;
2943: PetscFunctionReturn(PETSC_SUCCESS);
2944: }
2946: /*@C
2947: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
2948: processor's portion of the vector data. You MUST call `VecRestoreArray3d()`
2949: when you no longer need access to the array.
2951: Logically Collective
2953: Input Parameters:
2954: + x - the vector
2955: . m - first dimension of three dimensional array
2956: . n - second dimension of three dimensional array
2957: . p - third dimension of three dimensional array
2958: . mstart - first index you will use in first coordinate direction (often 0)
2959: . nstart - first index in the second coordinate direction (often 0)
2960: - pstart - first index in the third coordinate direction (often 0)
2962: Output Parameter:
2963: . a - location to put pointer to the array
2965: Level: developer
2967: Notes:
2968: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
2969: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2970: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2971: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
2973: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2975: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2976: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
2977: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2978: @*/
2979: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
2980: {
2981: PetscInt i, N, j;
2982: PetscScalar *aa, **b;
2984: PetscFunctionBegin;
2986: PetscAssertPointer(a, 8);
2988: PetscCall(VecGetLocalSize(x, &N));
2989: 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);
2990: PetscCall(VecGetArray(x, &aa));
2992: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
2993: b = (PetscScalar **)((*a) + m);
2994: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
2995: for (i = 0; i < m; i++)
2996: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset(aa, i * n * p + j * p - pstart);
2997: *a -= mstart;
2998: PetscFunctionReturn(PETSC_SUCCESS);
2999: }
3001: /*@C
3002: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
3003: processor's portion of the vector data. You MUST call `VecRestoreArray3dWrite()`
3004: when you no longer need access to the array.
3006: Logically Collective
3008: Input Parameters:
3009: + x - the vector
3010: . m - first dimension of three dimensional array
3011: . n - second dimension of three dimensional array
3012: . p - third dimension of three dimensional array
3013: . mstart - first index you will use in first coordinate direction (often 0)
3014: . nstart - first index in the second coordinate direction (often 0)
3015: - pstart - first index in the third coordinate direction (often 0)
3017: Output Parameter:
3018: . a - location to put pointer to the array
3020: Level: developer
3022: Notes:
3023: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3024: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3025: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3026: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3028: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3030: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3031: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3032: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3033: @*/
3034: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3035: {
3036: PetscInt i, N, j;
3037: PetscScalar *aa, **b;
3039: PetscFunctionBegin;
3041: PetscAssertPointer(a, 8);
3043: PetscCall(VecGetLocalSize(x, &N));
3044: 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);
3045: PetscCall(VecGetArrayWrite(x, &aa));
3047: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3048: b = (PetscScalar **)((*a) + m);
3049: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3050: for (i = 0; i < m; i++)
3051: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3053: *a -= mstart;
3054: PetscFunctionReturn(PETSC_SUCCESS);
3055: }
3057: /*@C
3058: VecRestoreArray3d - Restores a vector after `VecGetArray3d()` has been called.
3060: Logically Collective
3062: Input Parameters:
3063: + x - the vector
3064: . m - first dimension of three dimensional array
3065: . n - second dimension of the three dimensional array
3066: . p - third dimension of the three dimensional array
3067: . mstart - first index you will use in first coordinate direction (often 0)
3068: . nstart - first index in the second coordinate direction (often 0)
3069: . pstart - first index in the third coordinate direction (often 0)
3070: - a - location of pointer to array obtained from VecGetArray3d()
3072: Level: developer
3074: Notes:
3075: For regular PETSc vectors this routine does not involve any copies. For
3076: any special vectors that do not store local vector data in a contiguous
3077: array, this routine will copy the data back into the underlying
3078: vector data structure from the array obtained with `VecGetArray()`.
3080: This routine actually zeros out the `a` pointer.
3082: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3083: `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3084: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3085: @*/
3086: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3087: {
3088: void *dummy;
3090: PetscFunctionBegin;
3092: PetscAssertPointer(a, 8);
3094: dummy = (void *)(*a + mstart);
3095: PetscCall(PetscFree(dummy));
3096: PetscCall(VecRestoreArray(x, NULL));
3097: *a = NULL;
3098: PetscFunctionReturn(PETSC_SUCCESS);
3099: }
3101: /*@C
3102: VecRestoreArray3dWrite - Restores a vector after `VecGetArray3dWrite()` has been called.
3104: Logically Collective
3106: Input Parameters:
3107: + x - the vector
3108: . m - first dimension of three dimensional array
3109: . n - second dimension of the three dimensional array
3110: . p - third dimension of the three dimensional array
3111: . mstart - first index you will use in first coordinate direction (often 0)
3112: . nstart - first index in the second coordinate direction (often 0)
3113: . pstart - first index in the third coordinate direction (often 0)
3114: - a - location of pointer to array obtained from VecGetArray3d()
3116: Level: developer
3118: Notes:
3119: For regular PETSc vectors this routine does not involve any copies. For
3120: any special vectors that do not store local vector data in a contiguous
3121: array, this routine will copy the data back into the underlying
3122: vector data structure from the array obtained with `VecGetArray()`.
3124: This routine actually zeros out the `a` pointer.
3126: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3127: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3128: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3129: @*/
3130: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3131: {
3132: void *dummy;
3134: PetscFunctionBegin;
3136: PetscAssertPointer(a, 8);
3138: dummy = (void *)(*a + mstart);
3139: PetscCall(PetscFree(dummy));
3140: PetscCall(VecRestoreArrayWrite(x, NULL));
3141: *a = NULL;
3142: PetscFunctionReturn(PETSC_SUCCESS);
3143: }
3145: /*@C
3146: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3147: processor's portion of the vector data. You MUST call `VecRestoreArray4d()`
3148: when you no longer need access to the array.
3150: Logically Collective
3152: Input Parameters:
3153: + x - the vector
3154: . m - first dimension of four dimensional array
3155: . n - second dimension of four dimensional array
3156: . p - third dimension of four dimensional array
3157: . q - fourth dimension of four dimensional array
3158: . mstart - first index you will use in first coordinate direction (often 0)
3159: . nstart - first index in the second coordinate direction (often 0)
3160: . pstart - first index in the third coordinate direction (often 0)
3161: - qstart - first index in the fourth coordinate direction (often 0)
3163: Output Parameter:
3164: . a - location to put pointer to the array
3166: Level: developer
3168: Notes:
3169: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3170: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3171: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3172: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3174: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3176: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3177: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3178: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3179: @*/
3180: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3181: {
3182: PetscInt i, N, j, k;
3183: PetscScalar *aa, ***b, **c;
3185: PetscFunctionBegin;
3187: PetscAssertPointer(a, 10);
3189: PetscCall(VecGetLocalSize(x, &N));
3190: 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);
3191: PetscCall(VecGetArray(x, &aa));
3193: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3194: b = (PetscScalar ***)((*a) + m);
3195: c = (PetscScalar **)(b + m * n);
3196: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3197: for (i = 0; i < m; i++)
3198: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3199: for (i = 0; i < m; i++)
3200: for (j = 0; j < n; j++)
3201: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3202: *a -= mstart;
3203: PetscFunctionReturn(PETSC_SUCCESS);
3204: }
3206: /*@C
3207: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3208: processor's portion of the vector data. You MUST call `VecRestoreArray4dWrite()`
3209: when you no longer need access to the array.
3211: Logically Collective
3213: Input Parameters:
3214: + x - the vector
3215: . m - first dimension of four dimensional array
3216: . n - second dimension of four dimensional array
3217: . p - third dimension of four dimensional array
3218: . q - fourth dimension of four dimensional array
3219: . mstart - first index you will use in first coordinate direction (often 0)
3220: . nstart - first index in the second coordinate direction (often 0)
3221: . pstart - first index in the third coordinate direction (often 0)
3222: - qstart - first index in the fourth coordinate direction (often 0)
3224: Output Parameter:
3225: . a - location to put pointer to the array
3227: Level: developer
3229: Notes:
3230: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3231: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3232: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3233: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3235: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3237: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3238: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3239: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3240: @*/
3241: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3242: {
3243: PetscInt i, N, j, k;
3244: PetscScalar *aa, ***b, **c;
3246: PetscFunctionBegin;
3248: PetscAssertPointer(a, 10);
3250: PetscCall(VecGetLocalSize(x, &N));
3251: 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);
3252: PetscCall(VecGetArrayWrite(x, &aa));
3254: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3255: b = (PetscScalar ***)((*a) + m);
3256: c = (PetscScalar **)(b + m * n);
3257: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3258: for (i = 0; i < m; i++)
3259: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3260: for (i = 0; i < m; i++)
3261: for (j = 0; j < n; j++)
3262: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3263: *a -= mstart;
3264: PetscFunctionReturn(PETSC_SUCCESS);
3265: }
3267: /*@C
3268: VecRestoreArray4d - Restores a vector after `VecGetArray4d()` has been called.
3270: Logically Collective
3272: Input Parameters:
3273: + x - the vector
3274: . m - first dimension of four dimensional array
3275: . n - second dimension of the four dimensional array
3276: . p - third dimension of the four dimensional array
3277: . q - fourth dimension of the four dimensional array
3278: . mstart - first index you will use in first coordinate direction (often 0)
3279: . nstart - first index in the second coordinate direction (often 0)
3280: . pstart - first index in the third coordinate direction (often 0)
3281: . qstart - first index in the fourth coordinate direction (often 0)
3282: - a - location of pointer to array obtained from VecGetArray4d()
3284: Level: developer
3286: Notes:
3287: For regular PETSc vectors this routine does not involve any copies. For
3288: any special vectors that do not store local vector data in a contiguous
3289: array, this routine will copy the data back into the underlying
3290: vector data structure from the array obtained with `VecGetArray()`.
3292: This routine actually zeros out the `a` pointer.
3294: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3295: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3296: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`
3297: @*/
3298: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3299: {
3300: void *dummy;
3302: PetscFunctionBegin;
3304: PetscAssertPointer(a, 10);
3306: dummy = (void *)(*a + mstart);
3307: PetscCall(PetscFree(dummy));
3308: PetscCall(VecRestoreArray(x, NULL));
3309: *a = NULL;
3310: PetscFunctionReturn(PETSC_SUCCESS);
3311: }
3313: /*@C
3314: VecRestoreArray4dWrite - Restores a vector after `VecGetArray4dWrite()` has been called.
3316: Logically Collective
3318: Input Parameters:
3319: + x - the vector
3320: . m - first dimension of four dimensional array
3321: . n - second dimension of the four dimensional array
3322: . p - third dimension of the four dimensional array
3323: . q - fourth dimension of the four dimensional array
3324: . mstart - first index you will use in first coordinate direction (often 0)
3325: . nstart - first index in the second coordinate direction (often 0)
3326: . pstart - first index in the third coordinate direction (often 0)
3327: . qstart - first index in the fourth coordinate direction (often 0)
3328: - a - location of pointer to array obtained from `VecGetArray4d()`
3330: Level: developer
3332: Notes:
3333: For regular PETSc vectors this routine does not involve any copies. For
3334: any special vectors that do not store local vector data in a contiguous
3335: array, this routine will copy the data back into the underlying
3336: vector data structure from the array obtained with `VecGetArray()`.
3338: This routine actually zeros out the `a` pointer.
3340: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3341: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3342: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3343: @*/
3344: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3345: {
3346: void *dummy;
3348: PetscFunctionBegin;
3350: PetscAssertPointer(a, 10);
3352: dummy = (void *)(*a + mstart);
3353: PetscCall(PetscFree(dummy));
3354: PetscCall(VecRestoreArrayWrite(x, NULL));
3355: *a = NULL;
3356: PetscFunctionReturn(PETSC_SUCCESS);
3357: }
3359: /*@C
3360: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3361: processor's portion of the vector data. You MUST call `VecRestoreArray2dRead()`
3362: when you no longer need access to the array.
3364: Logically Collective
3366: Input Parameters:
3367: + x - the vector
3368: . m - first dimension of two dimensional array
3369: . n - second dimension of two dimensional array
3370: . mstart - first index you will use in first coordinate direction (often 0)
3371: - nstart - first index in the second coordinate direction (often 0)
3373: Output Parameter:
3374: . a - location to put pointer to the array
3376: Level: developer
3378: Notes:
3379: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
3380: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3381: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3382: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
3384: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3386: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3387: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3388: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3389: @*/
3390: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3391: {
3392: PetscInt i, N;
3393: const PetscScalar *aa;
3395: PetscFunctionBegin;
3397: PetscAssertPointer(a, 6);
3399: PetscCall(VecGetLocalSize(x, &N));
3400: 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);
3401: PetscCall(VecGetArrayRead(x, &aa));
3403: PetscCall(PetscMalloc1(m, a));
3404: for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3405: *a -= mstart;
3406: PetscFunctionReturn(PETSC_SUCCESS);
3407: }
3409: /*@C
3410: VecRestoreArray2dRead - Restores a vector after `VecGetArray2dRead()` has been called.
3412: Logically Collective
3414: Input Parameters:
3415: + x - the vector
3416: . m - first dimension of two dimensional array
3417: . n - second dimension of the two dimensional array
3418: . mstart - first index you will use in first coordinate direction (often 0)
3419: . nstart - first index in the second coordinate direction (often 0)
3420: - a - location of pointer to array obtained from VecGetArray2d()
3422: Level: developer
3424: Notes:
3425: For regular PETSc vectors this routine does not involve any copies. For
3426: any special vectors that do not store local vector data in a contiguous
3427: array, this routine will copy the data back into the underlying
3428: vector data structure from the array obtained with `VecGetArray()`.
3430: This routine actually zeros out the `a` pointer.
3432: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3433: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3434: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3435: @*/
3436: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3437: {
3438: void *dummy;
3440: PetscFunctionBegin;
3442: PetscAssertPointer(a, 6);
3444: dummy = (void *)(*a + mstart);
3445: PetscCall(PetscFree(dummy));
3446: PetscCall(VecRestoreArrayRead(x, NULL));
3447: *a = NULL;
3448: PetscFunctionReturn(PETSC_SUCCESS);
3449: }
3451: /*@C
3452: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3453: processor's portion of the vector data. You MUST call `VecRestoreArray1dRead()`
3454: when you no longer need access to the array.
3456: Logically Collective
3458: Input Parameters:
3459: + x - the vector
3460: . m - first dimension of two dimensional array
3461: - mstart - first index you will use in first coordinate direction (often 0)
3463: Output Parameter:
3464: . a - location to put pointer to the array
3466: Level: developer
3468: Notes:
3469: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3470: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3471: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3473: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3475: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3476: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3477: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3478: @*/
3479: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3480: {
3481: PetscInt N;
3483: PetscFunctionBegin;
3485: PetscAssertPointer(a, 4);
3487: PetscCall(VecGetLocalSize(x, &N));
3488: 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);
3489: PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3490: *a -= mstart;
3491: PetscFunctionReturn(PETSC_SUCCESS);
3492: }
3494: /*@C
3495: VecRestoreArray1dRead - Restores a vector after `VecGetArray1dRead()` has been called.
3497: Logically Collective
3499: Input Parameters:
3500: + x - the vector
3501: . m - first dimension of two dimensional array
3502: . mstart - first index you will use in first coordinate direction (often 0)
3503: - a - location of pointer to array obtained from `VecGetArray1dRead()`
3505: Level: developer
3507: Notes:
3508: For regular PETSc vectors this routine does not involve any copies. For
3509: any special vectors that do not store local vector data in a contiguous
3510: array, this routine will copy the data back into the underlying
3511: vector data structure from the array obtained with `VecGetArray1dRead()`.
3513: This routine actually zeros out the `a` pointer.
3515: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3516: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3517: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3518: @*/
3519: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3520: {
3521: PetscFunctionBegin;
3524: PetscCall(VecRestoreArrayRead(x, NULL));
3525: *a = NULL;
3526: PetscFunctionReturn(PETSC_SUCCESS);
3527: }
3529: /*@C
3530: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3531: processor's portion of the vector data. You MUST call `VecRestoreArray3dRead()`
3532: when you no longer need access to the array.
3534: Logically Collective
3536: Input Parameters:
3537: + x - the vector
3538: . m - first dimension of three dimensional array
3539: . n - second dimension of three dimensional array
3540: . p - third dimension of three dimensional array
3541: . mstart - first index you will use in first coordinate direction (often 0)
3542: . nstart - first index in the second coordinate direction (often 0)
3543: - pstart - first index in the third coordinate direction (often 0)
3545: Output Parameter:
3546: . a - location to put pointer to the array
3548: Level: developer
3550: Notes:
3551: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3552: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3553: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3554: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3dRead()`.
3556: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3558: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3559: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3560: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3561: @*/
3562: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3563: {
3564: PetscInt i, N, j;
3565: const PetscScalar *aa;
3566: PetscScalar **b;
3568: PetscFunctionBegin;
3570: PetscAssertPointer(a, 8);
3572: PetscCall(VecGetLocalSize(x, &N));
3573: 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);
3574: PetscCall(VecGetArrayRead(x, &aa));
3576: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3577: b = (PetscScalar **)((*a) + m);
3578: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3579: for (i = 0; i < m; i++)
3580: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset((PetscScalar *)aa, i * n * p + j * p - pstart);
3581: *a -= mstart;
3582: PetscFunctionReturn(PETSC_SUCCESS);
3583: }
3585: /*@C
3586: VecRestoreArray3dRead - Restores a vector after `VecGetArray3dRead()` has been called.
3588: Logically Collective
3590: Input Parameters:
3591: + x - the vector
3592: . m - first dimension of three dimensional array
3593: . n - second dimension of the three dimensional array
3594: . p - third dimension of the three dimensional array
3595: . mstart - first index you will use in first coordinate direction (often 0)
3596: . nstart - first index in the second coordinate direction (often 0)
3597: . pstart - first index in the third coordinate direction (often 0)
3598: - a - location of pointer to array obtained from `VecGetArray3dRead()`
3600: Level: developer
3602: Notes:
3603: For regular PETSc vectors this routine does not involve any copies. For
3604: any special vectors that do not store local vector data in a contiguous
3605: array, this routine will copy the data back into the underlying
3606: vector data structure from the array obtained with `VecGetArray()`.
3608: This routine actually zeros out the `a` pointer.
3610: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3611: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3612: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3613: @*/
3614: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3615: {
3616: void *dummy;
3618: PetscFunctionBegin;
3620: PetscAssertPointer(a, 8);
3622: dummy = (void *)(*a + mstart);
3623: PetscCall(PetscFree(dummy));
3624: PetscCall(VecRestoreArrayRead(x, NULL));
3625: *a = NULL;
3626: PetscFunctionReturn(PETSC_SUCCESS);
3627: }
3629: /*@C
3630: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3631: processor's portion of the vector data. You MUST call `VecRestoreArray4dRead()`
3632: when you no longer need access to the array.
3634: Logically Collective
3636: Input Parameters:
3637: + x - the vector
3638: . m - first dimension of four dimensional array
3639: . n - second dimension of four dimensional array
3640: . p - third dimension of four dimensional array
3641: . q - fourth dimension of four dimensional array
3642: . mstart - first index you will use in first coordinate direction (often 0)
3643: . nstart - first index in the second coordinate direction (often 0)
3644: . pstart - first index in the third coordinate direction (often 0)
3645: - qstart - first index in the fourth coordinate direction (often 0)
3647: Output Parameter:
3648: . a - location to put pointer to the array
3650: Level: beginner
3652: Notes:
3653: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3654: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3655: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3656: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3658: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3660: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3661: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3662: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3663: @*/
3664: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3665: {
3666: PetscInt i, N, j, k;
3667: const PetscScalar *aa;
3668: PetscScalar ***b, **c;
3670: PetscFunctionBegin;
3672: PetscAssertPointer(a, 10);
3674: PetscCall(VecGetLocalSize(x, &N));
3675: 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);
3676: PetscCall(VecGetArrayRead(x, &aa));
3678: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3679: b = (PetscScalar ***)((*a) + m);
3680: c = (PetscScalar **)(b + m * n);
3681: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3682: for (i = 0; i < m; i++)
3683: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3684: for (i = 0; i < m; i++)
3685: for (j = 0; j < n; j++)
3686: 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;
3687: *a -= mstart;
3688: PetscFunctionReturn(PETSC_SUCCESS);
3689: }
3691: /*@C
3692: VecRestoreArray4dRead - Restores a vector after `VecGetArray4d()` has been called.
3694: Logically Collective
3696: Input Parameters:
3697: + x - the vector
3698: . m - first dimension of four dimensional array
3699: . n - second dimension of the four dimensional array
3700: . p - third dimension of the four dimensional array
3701: . q - fourth dimension of the four dimensional array
3702: . mstart - first index you will use in first coordinate direction (often 0)
3703: . nstart - first index in the second coordinate direction (often 0)
3704: . pstart - first index in the third coordinate direction (often 0)
3705: . qstart - first index in the fourth coordinate direction (often 0)
3706: - a - location of pointer to array obtained from `VecGetArray4dRead()`
3708: Level: beginner
3710: Notes:
3711: For regular PETSc vectors this routine does not involve any copies. For
3712: any special vectors that do not store local vector data in a contiguous
3713: array, this routine will copy the data back into the underlying
3714: vector data structure from the array obtained with `VecGetArray()`.
3716: This routine actually zeros out the `a` pointer.
3718: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3719: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3720: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3721: @*/
3722: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3723: {
3724: void *dummy;
3726: PetscFunctionBegin;
3728: PetscAssertPointer(a, 10);
3730: dummy = (void *)(*a + mstart);
3731: PetscCall(PetscFree(dummy));
3732: PetscCall(VecRestoreArrayRead(x, NULL));
3733: *a = NULL;
3734: PetscFunctionReturn(PETSC_SUCCESS);
3735: }
3737: /*@
3738: VecLockGet - Get the current lock status of a vector
3740: Logically Collective
3742: Input Parameter:
3743: . x - the vector
3745: Output Parameter:
3746: . state - greater than zero indicates the vector is locked for read; less than zero indicates the vector is
3747: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3749: Level: advanced
3751: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3752: @*/
3753: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3754: {
3755: PetscFunctionBegin;
3757: PetscAssertPointer(state, 2);
3758: *state = x->lock;
3759: PetscFunctionReturn(PETSC_SUCCESS);
3760: }
3762: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3763: {
3764: PetscFunctionBegin;
3766: PetscAssertPointer(file, 2);
3767: PetscAssertPointer(func, 3);
3768: PetscAssertPointer(line, 4);
3769: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3770: {
3771: const int index = x->lockstack.currentsize - 1;
3773: *file = index < 0 ? NULL : x->lockstack.file[index];
3774: *func = index < 0 ? NULL : x->lockstack.function[index];
3775: *line = index < 0 ? 0 : x->lockstack.line[index];
3776: }
3777: #else
3778: *file = NULL;
3779: *func = NULL;
3780: *line = 0;
3781: #endif
3782: PetscFunctionReturn(PETSC_SUCCESS);
3783: }
3785: /*@
3786: VecLockReadPush - Push a read-only lock on a vector to prevent it from being written to
3788: Logically Collective
3790: Input Parameter:
3791: . x - the vector
3793: Level: intermediate
3795: Notes:
3796: If this is set then calls to `VecGetArray()` or `VecSetValues()` or any other routines that change the vectors values will generate an error.
3798: The call can be nested, i.e., called multiple times on the same vector, but each `VecLockReadPush()` has to have one matching
3799: `VecLockReadPop()`, which removes the latest read-only lock.
3801: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3802: @*/
3803: PetscErrorCode VecLockReadPush(Vec x)
3804: {
3805: PetscFunctionBegin;
3807: 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");
3808: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3809: {
3810: const char *file, *func;
3811: int index, line;
3813: if ((index = petscstack.currentsize - 2) < 0) {
3814: // vec was locked "outside" of petsc, either in user-land or main. the error message will
3815: // now show this function as the culprit, but it will include the stacktrace
3816: file = "unknown user-file";
3817: func = "unknown_user_function";
3818: line = 0;
3819: } else {
3820: file = petscstack.file[index];
3821: func = petscstack.function[index];
3822: line = petscstack.line[index];
3823: }
3824: PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
3825: }
3826: #endif
3827: PetscFunctionReturn(PETSC_SUCCESS);
3828: }
3830: /*@
3831: VecLockReadPop - Pop a read-only lock from a vector
3833: Logically Collective
3835: Input Parameter:
3836: . x - the vector
3838: Level: intermediate
3840: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
3841: @*/
3842: PetscErrorCode VecLockReadPop(Vec x)
3843: {
3844: PetscFunctionBegin;
3846: PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
3847: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3848: {
3849: const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];
3851: PetscStackPop_Private(x->lockstack, previous);
3852: }
3853: #endif
3854: PetscFunctionReturn(PETSC_SUCCESS);
3855: }
3857: /*@
3858: VecLockWriteSet - Lock or unlock a vector for exclusive read/write access
3860: Logically Collective
3862: Input Parameters:
3863: + x - the vector
3864: - flg - `PETSC_TRUE` to lock the vector for exclusive read/write access; `PETSC_FALSE` to unlock it.
3866: Level: intermediate
3868: Notes:
3869: The function is useful in split-phase computations, which usually have a begin phase and an end phase.
3870: One can call `VecLockWriteSet`(x,`PETSC_TRUE`) in the begin phase to lock a vector for exclusive
3871: access, and call `VecLockWriteSet`(x,`PETSC_FALSE`) in the end phase to unlock the vector from exclusive
3872: access. In this way, one is ensured no other operations can access the vector in between. The code may like
3874: .vb
3875: VecGetArray(x,&xdata); // begin phase
3876: VecLockWriteSet(v,PETSC_TRUE);
3878: Other operations, which can not access x anymore (they can access xdata, of course)
3880: VecRestoreArray(x,&vdata); // end phase
3881: VecLockWriteSet(v,PETSC_FALSE);
3882: .ve
3884: The call can not be nested on the same vector, in other words, one can not call `VecLockWriteSet`(x,`PETSC_TRUE`)
3885: again before calling `VecLockWriteSet`(v,`PETSC_FALSE`).
3887: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
3888: @*/
3889: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
3890: {
3891: PetscFunctionBegin;
3893: if (flg) {
3894: 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");
3895: 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");
3896: x->lock = -1;
3897: } else {
3898: 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");
3899: x->lock = 0;
3900: }
3901: PetscFunctionReturn(PETSC_SUCCESS);
3902: }