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/sfimpl.h"
6: #include "petscsystypes.h"
7: #include <petsc/private/vecimpl.h>
9: PetscInt VecGetSubVectorSavedStateId = -1;
11: #if PetscDefined(USE_DEBUG)
12: // this is a no-op '0' macro in optimized builds
13: PetscErrorCode VecValidValues_Internal(Vec vec, PetscInt argnum, PetscBool begin)
14: {
15: PetscFunctionBegin;
16: if (vec->petscnative || vec->ops->getarray) {
17: PetscInt n;
18: const PetscScalar *x;
19: PetscOffloadMask mask;
21: PetscCall(VecGetOffloadMask(vec, &mask));
22: if (!PetscOffloadHost(mask)) PetscFunctionReturn(PETSC_SUCCESS);
23: PetscCall(VecGetLocalSize(vec, &n));
24: PetscCall(VecGetArrayRead(vec, &x));
25: for (PetscInt i = 0; i < n; i++) {
26: if (begin) {
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 beginning of function: Parameter number %" PetscInt_FMT, i, argnum);
28: } else {
29: 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);
30: }
31: }
32: PetscCall(VecRestoreArrayRead(vec, &x));
33: }
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
36: #endif
38: /*@
39: VecMaxPointwiseDivide - Computes the maximum of the componentwise division `max = max_i abs(x[i]/y[i])`.
41: Logically Collective
43: Input Parameters:
44: + x - the numerators
45: - y - the denominators
47: Output Parameter:
48: . max - the result
50: Level: advanced
52: Notes:
53: `x` and `y` may be the same vector
55: if a particular `y[i]` is zero, it is treated as 1 in the above formula
57: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`
58: @*/
59: PetscErrorCode VecMaxPointwiseDivide(Vec x, Vec y, PetscReal *max)
60: {
61: PetscFunctionBegin;
64: PetscAssertPointer(max, 3);
67: PetscCheckSameTypeAndComm(x, 1, y, 2);
68: VecCheckSameSize(x, 1, y, 2);
69: VecCheckAssembled(x);
70: VecCheckAssembled(y);
71: PetscCall(VecLockReadPush(x));
72: PetscCall(VecLockReadPush(y));
73: PetscUseTypeMethod(x, maxpointwisedivide, y, max);
74: PetscCall(VecLockReadPop(x));
75: PetscCall(VecLockReadPop(y));
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: /*@
80: VecDot - Computes the vector dot product.
82: Collective
84: Input Parameters:
85: + x - first vector
86: - y - second vector
88: Output Parameter:
89: . val - the dot product
91: Level: intermediate
93: Notes for Users of Complex Numbers:
94: For complex vectors, `VecDot()` computes
95: $ val = (x,y) = y^H x,
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: $ val = (x,y) = y^T x,
102: where y^T denotes the transpose of y.
104: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDotRealPart()`
105: @*/
106: PetscErrorCode VecDot(Vec x, Vec y, PetscScalar *val)
107: {
108: PetscFunctionBegin;
111: PetscAssertPointer(val, 3);
114: PetscCheckSameTypeAndComm(x, 1, y, 2);
115: VecCheckSameSize(x, 1, y, 2);
116: VecCheckAssembled(x);
117: VecCheckAssembled(y);
119: PetscCall(VecLockReadPush(x));
120: PetscCall(VecLockReadPush(y));
121: PetscCall(PetscLogEventBegin(VEC_Dot, x, y, 0, 0));
122: PetscUseTypeMethod(x, dot, y, val);
123: PetscCall(PetscLogEventEnd(VEC_Dot, x, y, 0, 0));
124: PetscCall(VecLockReadPop(x));
125: PetscCall(VecLockReadPop(y));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: /*@
130: VecDotRealPart - Computes the real part of the vector dot product.
132: Collective
134: Input Parameters:
135: + x - first vector
136: - y - second vector
138: Output Parameter:
139: . val - the real part of the dot product;
141: Level: intermediate
143: Notes for Users of Complex Numbers:
144: See `VecDot()` for more details on the definition of the dot product for complex numbers
146: For real numbers this returns the same value as `VecDot()`
148: For complex numbers in C^n (that is a vector of n components with a complex number for each component) this is equal to the usual real dot product on the
149: the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)
151: Developer Notes:
152: This is not currently optimized to compute only the real part of the dot product.
154: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDot()`, `VecDotNorm2()`
155: @*/
156: PetscErrorCode VecDotRealPart(Vec x, Vec y, PetscReal *val)
157: {
158: PetscScalar fdot;
160: PetscFunctionBegin;
161: PetscCall(VecDot(x, y, &fdot));
162: *val = PetscRealPart(fdot);
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: /*@
167: VecNorm - Computes the vector norm.
169: Collective
171: Input Parameters:
172: + x - the vector
173: - type - the type of the norm requested
175: Output Parameter:
176: . val - the norm
178: Level: intermediate
180: Notes:
181: See `NormType` for descriptions of each norm.
183: For complex numbers `NORM_1` will return the traditional 1 norm of the 2 norm of the complex
184: numbers; that is the 1 norm of the absolute values of the complex entries. In PETSc 3.6 and
185: earlier releases it returned the 1 norm of the 1 norm of the complex entries (what is
186: returned by the BLAS routine `asum()`). Both are valid norms but most people expect the former.
188: This routine stashes the computed norm value, repeated calls before the vector entries are
189: changed are then rapid since the precomputed value is immediately available. Certain vector
190: operations such as `VecSet()` store the norms so the value is immediately available and does
191: not need to be explicitly computed. `VecScale()` updates any stashed norm values, thus calls
192: after `VecScale()` do not need to explicitly recompute the norm.
194: .seealso: [](ch_vectors), `Vec`, `NormType`, `VecDot()`, `VecTDot()`, `VecDotBegin()`, `VecDotEnd()`, `VecNormAvailable()`,
195: `VecNormBegin()`, `VecNormEnd()`, `NormType()`
196: @*/
197: PetscErrorCode VecNorm(Vec x, NormType type, PetscReal *val)
198: {
199: PetscBool flg = PETSC_TRUE;
201: PetscFunctionBegin;
204: VecCheckAssembled(x);
206: PetscAssertPointer(val, 3);
208: PetscCall(VecNormAvailable(x, type, &flg, val));
209: // check that all MPI processes call this routine together and have same availability
210: if (PetscDefined(USE_DEBUG)) {
211: PetscBool minflg;
213: PetscCall(MPIU_Allreduce(&flg, &minflg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)x)));
214: PetscCheck(flg == minflg, PetscObjectComm((PetscObject)(x)), PETSC_ERR_ARG_WRONGSTATE, "Some MPI processes have cached norm, others do not. This may happen when some MPI processes call VecGetArray() and some others do not.");
215: }
216: if (flg) PetscFunctionReturn(PETSC_SUCCESS);
218: PetscCall(VecLockReadPush(x));
219: PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
220: PetscUseTypeMethod(x, norm, type, val);
221: PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
222: PetscCall(VecLockReadPop(x));
224: if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: /*@
229: VecNormAvailable - Returns the vector norm if it is already known. That is, it has been previously computed and cached in the vector
231: Not Collective
233: Input Parameters:
234: + x - the vector
235: - 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
236: `NORM_1_AND_2`, which computes both norms and stores them
237: in a two element array.
239: Output Parameters:
240: + available - `PETSC_TRUE` if the val returned is valid
241: - val - the norm
243: Level: intermediate
245: Developer Notes:
246: `PETSC_HAVE_SLOW_BLAS_NORM2` will cause a C (loop unrolled) version of the norm to be used, rather
247: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
248: (as opposed to vendor provided) because the FORTRAN BLAS `NRM2()` routine is very slow.
250: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
251: `VecNormBegin()`, `VecNormEnd()`
252: @*/
253: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
254: {
255: PetscFunctionBegin;
258: PetscAssertPointer(available, 3);
259: PetscAssertPointer(val, 4);
261: if (type == NORM_1_AND_2) {
262: *available = PETSC_FALSE;
263: } else {
264: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
265: }
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: /*@
270: VecNormalize - Normalizes a vector by its 2-norm.
272: Collective
274: Input Parameter:
275: . x - the vector
277: Output Parameter:
278: . val - the vector norm before normalization. May be `NULL` if the value is not needed.
280: Level: intermediate
282: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
283: @*/
284: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
285: {
286: PetscReal norm;
288: PetscFunctionBegin;
291: PetscCall(VecSetErrorIfLocked(x, 1));
292: if (val) PetscAssertPointer(val, 2);
293: PetscCall(PetscLogEventBegin(VEC_Normalize, x, 0, 0, 0));
294: PetscCall(VecNorm(x, NORM_2, &norm));
295: if (norm == 0.0) PetscCall(PetscInfo(x, "Vector of zero norm can not be normalized; Returning only the zero norm\n"));
296: else if (PetscIsInfOrNanReal(norm)) PetscCall(PetscInfo(x, "Vector with Inf or Nan norm can not be normalized; Returning only the norm\n"));
297: else {
298: PetscScalar s = 1.0 / norm;
299: PetscCall(VecScale(x, s));
300: }
301: PetscCall(PetscLogEventEnd(VEC_Normalize, x, 0, 0, 0));
302: if (val) *val = norm;
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: /*@C
307: VecMax - Determines the vector component with maximum real part and its location.
309: Collective
311: Input Parameter:
312: . x - the vector
314: Output Parameters:
315: + p - the index of `val` (pass `NULL` if you don't want this) in the vector
316: - val - the maximum component
318: Level: intermediate
320: Notes:
321: Returns the value `PETSC_MIN_REAL` and negative `p` if the vector is of length 0.
323: Returns the smallest index with the maximum value
325: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
326: @*/
327: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
328: {
329: PetscFunctionBegin;
332: VecCheckAssembled(x);
333: if (p) PetscAssertPointer(p, 2);
334: PetscAssertPointer(val, 3);
335: PetscCall(VecLockReadPush(x));
336: PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
337: PetscUseTypeMethod(x, max, p, val);
338: PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
339: PetscCall(VecLockReadPop(x));
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: /*@C
344: VecMin - Determines the vector component with minimum real part and its location.
346: Collective
348: Input Parameter:
349: . x - the vector
351: Output Parameters:
352: + p - the index of `val` (pass `NULL` if you don't want this location) in the vector
353: - val - the minimum component
355: Level: intermediate
357: Notes:
358: Returns the value `PETSC_MAX_REAL` and negative `p` if the vector is of length 0.
360: This returns the smallest index with the minimum value
362: .seealso: [](ch_vectors), `Vec`, `VecMax()`
363: @*/
364: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
365: {
366: PetscFunctionBegin;
369: VecCheckAssembled(x);
370: if (p) PetscAssertPointer(p, 2);
371: PetscAssertPointer(val, 3);
372: PetscCall(VecLockReadPush(x));
373: PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
374: PetscUseTypeMethod(x, min, p, val);
375: PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
376: PetscCall(VecLockReadPop(x));
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: /*@
381: VecTDot - Computes an indefinite vector dot product. That is, this
382: routine does NOT use the complex conjugate.
384: Collective
386: Input Parameters:
387: + x - first vector
388: - y - second vector
390: Output Parameter:
391: . val - the dot product
393: Level: intermediate
395: Notes for Users of Complex Numbers:
396: For complex vectors, `VecTDot()` computes the indefinite form
397: $ val = (x,y) = y^T x,
398: where y^T denotes the transpose of y.
400: Use `VecDot()` for the inner product
401: $ val = (x,y) = y^H x,
402: where y^H denotes the conjugate transpose of y.
404: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
405: @*/
406: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
407: {
408: PetscFunctionBegin;
411: PetscAssertPointer(val, 3);
414: PetscCheckSameTypeAndComm(x, 1, y, 2);
415: VecCheckSameSize(x, 1, y, 2);
416: VecCheckAssembled(x);
417: VecCheckAssembled(y);
419: PetscCall(VecLockReadPush(x));
420: PetscCall(VecLockReadPush(y));
421: PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
422: PetscUseTypeMethod(x, tdot, y, val);
423: PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
424: PetscCall(VecLockReadPop(x));
425: PetscCall(VecLockReadPop(y));
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: PetscErrorCode VecScaleAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
430: {
431: PetscReal norms[4];
432: PetscBool flgs[4];
433: PetscScalar one = 1.0;
435: PetscFunctionBegin;
438: VecCheckAssembled(x);
439: PetscCall(VecSetErrorIfLocked(x, 1));
440: if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);
442: /* get current stashed norms */
443: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
445: PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
446: VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
447: PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));
449: PetscCall(PetscObjectStateIncrease((PetscObject)x));
450: /* put the scaled stashed norms back into the Vec */
451: for (PetscInt i = 0; i < 4; i++) {
452: PetscReal ar = PetscAbsScalar(alpha);
453: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
454: }
455: PetscFunctionReturn(PETSC_SUCCESS);
456: }
458: /*@
459: VecScale - Scales a vector.
461: Not Collective
463: Input Parameters:
464: + x - the vector
465: - alpha - the scalar
467: Level: intermediate
469: Note:
470: For a vector with n components, `VecScale()` computes x[i] = alpha * x[i], for i=1,...,n.
472: .seealso: [](ch_vectors), `Vec`, `VecSet()`
473: @*/
474: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
475: {
476: PetscFunctionBegin;
477: PetscCall(VecScaleAsync_Private(x, alpha, NULL));
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
482: {
483: PetscFunctionBegin;
486: VecCheckAssembled(x);
488: PetscCall(VecSetErrorIfLocked(x, 1));
490: if (alpha == 0) {
491: PetscReal norm;
492: PetscBool set;
494: PetscCall(VecNormAvailable(x, NORM_2, &set, &norm));
495: if (set == PETSC_TRUE && norm == 0) PetscFunctionReturn(PETSC_SUCCESS);
496: }
497: PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
498: VecMethodDispatch(x, dctx, VecAsyncFnName(Set), set, (Vec, PetscScalar, PetscDeviceContext), alpha);
499: PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
500: PetscCall(PetscObjectStateIncrease((PetscObject)x));
502: /* norms can be simply set (if |alpha|*N not too large) */
504: {
505: PetscReal val = PetscAbsScalar(alpha);
506: const PetscInt N = x->map->N;
508: if (N == 0) {
509: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0l));
510: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
511: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
512: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
513: } else if (val > PETSC_MAX_REAL / N) {
514: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
515: } else {
516: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
517: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
518: val *= PetscSqrtReal((PetscReal)N);
519: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
520: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
521: }
522: }
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: /*@
527: VecSet - Sets all components of a vector to a single scalar value.
529: Logically Collective
531: Input Parameters:
532: + x - the vector
533: - alpha - the scalar
535: Level: beginner
537: Notes:
538: For a vector of dimension n, `VecSet()` sets x[i] = alpha, for i=1,...,n,
539: so that all vector entries then equal the identical
540: scalar value, `alpha`. Use the more general routine
541: `VecSetValues()` to set different vector entries.
543: You CANNOT call this after you have called `VecSetValues()` but before you call
544: `VecAssemblyBegin()`
546: If `alpha` is zero and the norm of the vector is known to be zero then this skips the unneeded zeroing process
548: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
549: @*/
550: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
551: {
552: PetscFunctionBegin;
553: PetscCall(VecSetAsync_Private(x, alpha, NULL));
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }
557: PetscErrorCode VecAXPYAsync_Private(Vec y, PetscScalar alpha, Vec x, PetscDeviceContext dctx)
558: {
559: PetscFunctionBegin;
564: PetscCheckSameTypeAndComm(x, 3, y, 1);
565: VecCheckSameSize(x, 3, y, 1);
566: VecCheckAssembled(x);
567: VecCheckAssembled(y);
569: if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
570: PetscCall(VecSetErrorIfLocked(y, 1));
571: if (x == y) {
572: PetscCall(VecScale(y, alpha + 1.0));
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
575: PetscCall(VecLockReadPush(x));
576: PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
577: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPY), axpy, (Vec, PetscScalar, Vec, PetscDeviceContext), alpha, x);
578: PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
579: PetscCall(VecLockReadPop(x));
580: PetscCall(PetscObjectStateIncrease((PetscObject)y));
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
583: /*@
584: VecAXPY - Computes `y = alpha x + y`.
586: Logically Collective
588: Input Parameters:
589: + alpha - the scalar
590: . x - vector scale by `alpha`
591: - y - vector accumulated into
593: Output Parameter:
594: . y - output vector
596: Level: intermediate
598: Notes:
599: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
600: .vb
601: VecAXPY(y,alpha,x) y = alpha x + y
602: VecAYPX(y,beta,x) y = x + beta y
603: VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
604: VecWAXPY(w,alpha,x,y) w = alpha x + y
605: VecAXPBYPCZ(z,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
606: VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
607: .ve
609: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
610: @*/
611: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
612: {
613: PetscFunctionBegin;
614: PetscCall(VecAXPYAsync_Private(y, alpha, x, NULL));
615: PetscFunctionReturn(PETSC_SUCCESS);
616: }
618: PetscErrorCode VecAYPXAsync_Private(Vec y, PetscScalar beta, Vec x, PetscDeviceContext dctx)
619: {
620: PetscFunctionBegin;
625: PetscCheckSameTypeAndComm(x, 3, y, 1);
626: VecCheckSameSize(x, 1, y, 3);
627: VecCheckAssembled(x);
628: VecCheckAssembled(y);
630: PetscCall(VecSetErrorIfLocked(y, 1));
631: if (x == y) {
632: PetscCall(VecScale(y, beta + 1.0));
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
635: PetscCall(VecLockReadPush(x));
636: if (beta == (PetscScalar)0.0) {
637: PetscCall(VecCopy(x, y));
638: } else {
639: PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
640: VecMethodDispatch(y, dctx, VecAsyncFnName(AYPX), aypx, (Vec, PetscScalar, Vec, PetscDeviceContext), beta, x);
641: PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
642: PetscCall(PetscObjectStateIncrease((PetscObject)y));
643: }
644: PetscCall(VecLockReadPop(x));
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
648: /*@
649: VecAYPX - Computes `y = x + beta y`.
651: Logically Collective
653: Input Parameters:
654: + beta - the scalar
655: . x - the unscaled vector
656: - y - the vector to be scaled
658: Output Parameter:
659: . y - output vector
661: Level: intermediate
663: Developer Notes:
664: The implementation is optimized for `beta` of -1.0, 0.0, and 1.0
666: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
667: @*/
668: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
669: {
670: PetscFunctionBegin;
671: PetscCall(VecAYPXAsync_Private(y, beta, x, NULL));
672: PetscFunctionReturn(PETSC_SUCCESS);
673: }
675: PetscErrorCode VecAXPBYAsync_Private(Vec y, PetscScalar alpha, PetscScalar beta, Vec x, PetscDeviceContext dctx)
676: {
677: PetscFunctionBegin;
682: PetscCheckSameTypeAndComm(x, 4, y, 1);
683: VecCheckSameSize(y, 1, x, 4);
684: VecCheckAssembled(x);
685: VecCheckAssembled(y);
688: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
689: if (x == y) {
690: PetscCall(VecScale(y, alpha + beta));
691: PetscFunctionReturn(PETSC_SUCCESS);
692: }
694: PetscCall(VecSetErrorIfLocked(y, 1));
695: PetscCall(VecLockReadPush(x));
696: PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
697: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPBY), axpby, (Vec, PetscScalar, PetscScalar, Vec, PetscDeviceContext), alpha, beta, x);
698: PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
699: PetscCall(PetscObjectStateIncrease((PetscObject)y));
700: PetscCall(VecLockReadPop(x));
701: PetscFunctionReturn(PETSC_SUCCESS);
702: }
703: /*@
704: VecAXPBY - Computes `y = alpha x + beta y`.
706: Logically Collective
708: Input Parameters:
709: + alpha - first scalar
710: . beta - second scalar
711: . x - the first scaled vector
712: - y - the second scaled vector
714: Output Parameter:
715: . y - output vector
717: Level: intermediate
719: Developer Notes:
720: The implementation is optimized for `alpha` and/or `beta` values of 0.0 and 1.0
722: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
723: @*/
724: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
725: {
726: PetscFunctionBegin;
727: PetscCall(VecAXPBYAsync_Private(y, alpha, beta, x, NULL));
728: PetscFunctionReturn(PETSC_SUCCESS);
729: }
731: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
732: {
733: PetscFunctionBegin;
740: PetscCheckSameTypeAndComm(x, 5, y, 6);
741: PetscCheckSameTypeAndComm(x, 5, z, 1);
742: VecCheckSameSize(x, 1, y, 5);
743: VecCheckSameSize(x, 1, z, 6);
744: PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
745: PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
746: VecCheckAssembled(x);
747: VecCheckAssembled(y);
748: VecCheckAssembled(z);
752: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
754: PetscCall(VecSetErrorIfLocked(z, 1));
755: PetscCall(VecLockReadPush(x));
756: PetscCall(VecLockReadPush(y));
757: PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
758: VecMethodDispatch(z, dctx, VecAsyncFnName(AXPBYPCZ), axpbypcz, (Vec, PetscScalar, PetscScalar, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, beta, gamma, x, y);
759: PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
760: PetscCall(PetscObjectStateIncrease((PetscObject)z));
761: PetscCall(VecLockReadPop(x));
762: PetscCall(VecLockReadPop(y));
763: PetscFunctionReturn(PETSC_SUCCESS);
764: }
765: /*@
766: VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`
768: Logically Collective
770: Input Parameters:
771: + alpha - first scalar
772: . beta - second scalar
773: . gamma - third scalar
774: . x - first vector
775: . y - second vector
776: - z - third vector
778: Output Parameter:
779: . z - output vector
781: Level: intermediate
783: Note:
784: `x`, `y` and `z` must be different vectors
786: Developer Notes:
787: The implementation is optimized for `alpha` of 1.0 and `gamma` of 1.0 or 0.0
789: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
790: @*/
791: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
792: {
793: PetscFunctionBegin;
794: PetscCall(VecAXPBYPCZAsync_Private(z, alpha, beta, gamma, x, y, NULL));
795: PetscFunctionReturn(PETSC_SUCCESS);
796: }
798: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
799: {
800: PetscFunctionBegin;
807: PetscCheckSameTypeAndComm(x, 3, y, 4);
808: PetscCheckSameTypeAndComm(y, 4, w, 1);
809: VecCheckSameSize(x, 3, y, 4);
810: VecCheckSameSize(x, 3, w, 1);
811: PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
812: PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
813: VecCheckAssembled(x);
814: VecCheckAssembled(y);
816: PetscCall(VecSetErrorIfLocked(w, 1));
818: PetscCall(VecLockReadPush(x));
819: PetscCall(VecLockReadPush(y));
820: if (alpha == (PetscScalar)0.0) {
821: PetscCall(VecCopyAsync_Private(y, w, dctx));
822: } else {
823: PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
824: VecMethodDispatch(w, dctx, VecAsyncFnName(WAXPY), waxpy, (Vec, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, x, y);
825: PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
826: PetscCall(PetscObjectStateIncrease((PetscObject)w));
827: }
828: PetscCall(VecLockReadPop(x));
829: PetscCall(VecLockReadPop(y));
830: PetscFunctionReturn(PETSC_SUCCESS);
831: }
833: /*@
834: VecWAXPY - Computes `w = alpha x + y`.
836: Logically Collective
838: Input Parameters:
839: + alpha - the scalar
840: . x - first vector, multiplied by `alpha`
841: - y - second vector
843: Output Parameter:
844: . w - the result
846: Level: intermediate
848: Note:
849: `w` cannot be either `x` or `y`, but `x` and `y` can be the same
851: Developer Notes:
852: The implementation is optimized for alpha of -1.0, 0.0, and 1.0
854: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
855: @*/
856: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
857: {
858: PetscFunctionBegin;
859: PetscCall(VecWAXPYAsync_Private(w, alpha, x, y, NULL));
860: PetscFunctionReturn(PETSC_SUCCESS);
861: }
863: /*@C
864: VecSetValues - Inserts or adds values into certain locations of a vector.
866: Not Collective
868: Input Parameters:
869: + x - vector to insert in
870: . ni - number of elements to add
871: . ix - indices where to add
872: . y - array of values
873: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries
875: Level: beginner
877: Notes:
878: .vb
879: `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
880: .ve
882: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
883: options cannot be mixed without intervening calls to the assembly
884: routines.
886: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
887: MUST be called after all calls to `VecSetValues()` have been completed.
889: VecSetValues() uses 0-based indices in Fortran as well as in C.
891: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
892: negative indices may be passed in ix. These rows are
893: simply ignored. This allows easily inserting element load matrices
894: with homogeneous Dirichlet boundary conditions that you don't want represented
895: in the vector.
897: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
898: `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
899: @*/
900: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
901: {
902: PetscFunctionBeginHot;
904: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
905: PetscAssertPointer(ix, 3);
906: PetscAssertPointer(y, 4);
909: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
910: PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
911: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
912: PetscCall(PetscObjectStateIncrease((PetscObject)x));
913: PetscFunctionReturn(PETSC_SUCCESS);
914: }
916: /*@C
917: VecGetValues - Gets values from certain locations of a vector. Currently
918: can only get values on the same processor on which they are owned
920: Not Collective
922: Input Parameters:
923: + x - vector to get values from
924: . ni - number of elements to get
925: - ix - indices where to get them from (in global 1d numbering)
927: Output Parameter:
928: . y - array of values
930: Level: beginner
932: Notes:
933: The user provides the allocated array y; it is NOT allocated in this routine
935: `VecGetValues()` gets y[i] = x[ix[i]], for i=0,...,ni-1.
937: `VecAssemblyBegin()` and `VecAssemblyEnd()` MUST be called before calling this if `VecSetValues()` or related routine has been called
939: VecGetValues() uses 0-based indices in Fortran as well as in C.
941: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
942: negative indices may be passed in ix. These rows are
943: simply ignored.
945: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
946: @*/
947: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
948: {
949: PetscFunctionBegin;
951: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
952: PetscAssertPointer(ix, 3);
953: PetscAssertPointer(y, 4);
955: VecCheckAssembled(x);
956: PetscUseTypeMethod(x, getvalues, ni, ix, y);
957: PetscFunctionReturn(PETSC_SUCCESS);
958: }
960: /*@C
961: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
963: Not Collective
965: Input Parameters:
966: + x - vector to insert in
967: . ni - number of blocks to add
968: . ix - indices where to add in block count, rather than element count
969: . y - array of values
970: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries
972: Level: intermediate
974: Notes:
975: `VecSetValuesBlocked()` sets x[bs*ix[i]+j] = y[bs*i+j],
976: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
978: Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
979: options cannot be mixed without intervening calls to the assembly
980: routines.
982: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
983: MUST be called after all calls to `VecSetValuesBlocked()` have been completed.
985: `VecSetValuesBlocked()` uses 0-based indices in Fortran as well as in C.
987: Negative indices may be passed in ix, these rows are
988: simply ignored. This allows easily inserting element load matrices
989: with homogeneous Dirichlet boundary conditions that you don't want represented
990: in the vector.
992: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
993: `VecSetValues()`
994: @*/
995: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
996: {
997: PetscFunctionBeginHot;
999: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1000: PetscAssertPointer(ix, 3);
1001: PetscAssertPointer(y, 4);
1004: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1005: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1006: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1007: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1008: PetscFunctionReturn(PETSC_SUCCESS);
1009: }
1011: /*@C
1012: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1013: using a local ordering of the nodes.
1015: Not Collective
1017: Input Parameters:
1018: + x - vector to insert in
1019: . ni - number of elements to add
1020: . ix - indices where to add
1021: . y - array of values
1022: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1024: Level: intermediate
1026: Notes:
1027: `VecSetValuesLocal()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
1029: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
1030: options cannot be mixed without intervening calls to the assembly
1031: routines.
1033: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1034: MUST be called after all calls to `VecSetValuesLocal()` have been completed.
1036: `VecSetValuesLocal()` uses 0-based indices in Fortran as well as in C.
1038: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1039: `VecSetValuesBlockedLocal()`
1040: @*/
1041: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1042: {
1043: PetscInt lixp[128], *lix = lixp;
1045: PetscFunctionBeginHot;
1047: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1048: PetscAssertPointer(ix, 3);
1049: PetscAssertPointer(y, 4);
1052: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1053: if (!x->ops->setvalueslocal) {
1054: if (x->map->mapping) {
1055: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1056: PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1057: PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1058: if (ni > 128) PetscCall(PetscFree(lix));
1059: } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1060: } else PetscUseTypeMethod(x, setvalueslocal, ni, ix, y, iora);
1061: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1062: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1063: PetscFunctionReturn(PETSC_SUCCESS);
1064: }
1066: /*@
1067: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1068: using a local ordering of the nodes.
1070: Not Collective
1072: Input Parameters:
1073: + x - vector to insert in
1074: . ni - number of blocks to add
1075: . ix - indices where to add in block count, not element count
1076: . y - array of values
1077: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1079: Level: intermediate
1081: Notes:
1082: `VecSetValuesBlockedLocal()` sets x[bs*ix[i]+j] = y[bs*i+j],
1083: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with `VecSetBlockSize()`.
1085: Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1086: options cannot be mixed without intervening calls to the assembly
1087: routines.
1089: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1090: MUST be called after all calls to `VecSetValuesBlockedLocal()` have been completed.
1092: `VecSetValuesBlockedLocal()` uses 0-based indices in Fortran as well as in C.
1094: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1095: `VecSetLocalToGlobalMapping()`
1096: @*/
1097: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1098: {
1099: PetscInt lixp[128], *lix = lixp;
1101: PetscFunctionBeginHot;
1103: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1104: PetscAssertPointer(ix, 3);
1105: PetscAssertPointer(y, 4);
1107: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1108: if (x->map->mapping) {
1109: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1110: PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1111: PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1112: if (ni > 128) PetscCall(PetscFree(lix));
1113: } else {
1114: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1115: }
1116: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1117: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1118: PetscFunctionReturn(PETSC_SUCCESS);
1119: }
1121: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1122: {
1123: PetscFunctionBegin;
1126: VecCheckAssembled(x);
1128: if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1129: PetscAssertPointer(y, 3);
1130: for (PetscInt i = 0; i < nv; ++i) {
1133: PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1134: VecCheckSameSize(x, 1, y[i], 3);
1135: VecCheckAssembled(y[i]);
1136: PetscCall(VecLockReadPush(y[i]));
1137: }
1138: PetscAssertPointer(result, 4);
1141: PetscCall(VecLockReadPush(x));
1142: PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1143: PetscCall((*mxdot)(x, nv, y, result));
1144: PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1145: PetscCall(VecLockReadPop(x));
1146: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1147: PetscFunctionReturn(PETSC_SUCCESS);
1148: }
1150: /*@
1151: VecMTDot - Computes indefinite vector multiple dot products.
1152: That is, it does NOT use the complex conjugate.
1154: Collective
1156: Input Parameters:
1157: + x - one vector
1158: . nv - number of vectors
1159: - y - array of vectors. Note that vectors are pointers
1161: Output Parameter:
1162: . val - array of the dot products
1164: Level: intermediate
1166: Notes for Users of Complex Numbers:
1167: For complex vectors, `VecMTDot()` computes the indefinite form
1168: $ val = (x,y) = y^T x,
1169: where y^T denotes the transpose of y.
1171: Use `VecMDot()` for the inner product
1172: $ val = (x,y) = y^H x,
1173: where y^H denotes the conjugate transpose of y.
1175: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1176: @*/
1177: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1178: {
1179: PetscFunctionBegin;
1181: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1182: PetscFunctionReturn(PETSC_SUCCESS);
1183: }
1185: /*@
1186: VecMDot - Computes multiple vector dot products.
1188: Collective
1190: Input Parameters:
1191: + x - one vector
1192: . nv - number of vectors
1193: - y - array of vectors.
1195: Output Parameter:
1196: . val - array of the dot products (does not allocate the array)
1198: Level: intermediate
1200: Notes for Users of Complex Numbers:
1201: For complex vectors, `VecMDot()` computes
1202: $ val = (x,y) = y^H x,
1203: where y^H denotes the conjugate transpose of y.
1205: Use `VecMTDot()` for the indefinite form
1206: $ val = (x,y) = y^T x,
1207: where y^T denotes the transpose of y.
1209: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`
1210: @*/
1211: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1212: {
1213: PetscFunctionBegin;
1215: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1216: PetscFunctionReturn(PETSC_SUCCESS);
1217: }
1219: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
1220: {
1221: PetscFunctionBegin;
1223: VecCheckAssembled(y);
1225: PetscCall(VecSetErrorIfLocked(y, 1));
1226: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1227: if (nv) {
1228: PetscInt zeros = 0;
1230: PetscAssertPointer(alpha, 3);
1231: PetscAssertPointer(x, 4);
1232: for (PetscInt i = 0; i < nv; ++i) {
1236: PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1237: VecCheckSameSize(y, 1, x[i], 4);
1238: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1239: VecCheckAssembled(x[i]);
1240: PetscCall(VecLockReadPush(x[i]));
1241: zeros += alpha[i] == (PetscScalar)0.0;
1242: }
1244: if (zeros < nv) {
1245: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1246: VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1247: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1248: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1249: }
1251: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1252: }
1253: PetscFunctionReturn(PETSC_SUCCESS);
1254: }
1256: /*@
1257: VecMAXPY - Computes `y = y + sum alpha[i] x[i]`
1259: Logically Collective
1261: Input Parameters:
1262: + nv - number of scalars and x-vectors
1263: . alpha - array of scalars
1264: . y - one vector
1265: - x - array of vectors
1267: Level: intermediate
1269: Note:
1270: `y` cannot be any of the `x` vectors
1272: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`,`VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1273: @*/
1274: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1275: {
1276: PetscFunctionBegin;
1277: PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1278: PetscFunctionReturn(PETSC_SUCCESS);
1279: }
1281: /*@
1282: VecMAXPBY - Computes `y = beta y + sum alpha[i] x[i]`
1284: Logically Collective
1286: Input Parameters:
1287: + nv - number of scalars and x-vectors
1288: . alpha - array of scalars
1289: . beta - scalar
1290: . y - one vector
1291: - x - array of vectors
1293: Level: intermediate
1295: Note:
1296: `y` cannot be any of the `x` vectors.
1298: Developer Notes:
1299: This is a convenience routine, but implementations might be able to optimize it, for example, when `beta` is zero.
1301: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1302: @*/
1303: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
1304: {
1305: PetscFunctionBegin;
1307: VecCheckAssembled(y);
1309: PetscCall(VecSetErrorIfLocked(y, 1));
1310: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1313: if (y->ops->maxpby) {
1314: PetscInt zeros = 0;
1316: if (nv) {
1317: PetscAssertPointer(alpha, 3);
1318: PetscAssertPointer(x, 5);
1319: }
1321: for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
1325: PetscCheckSameTypeAndComm(y, 1, x[i], 5);
1326: VecCheckSameSize(y, 1, x[i], 5);
1327: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1328: VecCheckAssembled(x[i]);
1329: PetscCall(VecLockReadPush(x[i]));
1330: zeros += alpha[i] == (PetscScalar)0.0;
1331: }
1333: if (zeros < nv) { // has nonzero alpha
1334: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1335: PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1336: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1337: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1338: } else {
1339: PetscCall(VecScale(y, beta));
1340: }
1342: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1343: } else { // no maxpby
1344: if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1345: else PetscCall(VecScale(y, beta));
1346: PetscCall(VecMAXPY(y, nv, alpha, x));
1347: }
1348: PetscFunctionReturn(PETSC_SUCCESS);
1349: }
1351: /*@
1352: VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
1353: in the order they appear in the array. The concatenated vector resides on the same
1354: communicator and is the same type as the source vectors.
1356: Collective
1358: Input Parameters:
1359: + nx - number of vectors to be concatenated
1360: - X - array containing the vectors to be concatenated in the order of concatenation
1362: Output Parameters:
1363: + Y - concatenated vector
1364: - x_is - array of index sets corresponding to the concatenated components of `Y` (pass `NULL` if not needed)
1366: Level: advanced
1368: Notes:
1369: Concatenation is similar to the functionality of a `VECNEST` object; they both represent combination of
1370: different vector spaces. However, concatenated vectors do not store any information about their
1371: sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1372: manipulation of data in the concatenated vector that corresponds to the original components at creation.
1374: This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1375: has to operate on combined vector spaces and cannot utilize `VECNEST` objects due to incompatibility with
1376: bound projections.
1378: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1379: @*/
1380: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1381: {
1382: MPI_Comm comm;
1383: VecType vec_type;
1384: Vec Ytmp, Xtmp;
1385: IS *is_tmp;
1386: PetscInt i, shift = 0, Xnl, Xng, Xbegin;
1388: PetscFunctionBegin;
1392: PetscAssertPointer(Y, 3);
1394: if ((*X)->ops->concatenate) {
1395: /* use the dedicated concatenation function if available */
1396: PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1397: } else {
1398: /* loop over vectors and start creating IS */
1399: comm = PetscObjectComm((PetscObject)(*X));
1400: PetscCall(VecGetType(*X, &vec_type));
1401: PetscCall(PetscMalloc1(nx, &is_tmp));
1402: for (i = 0; i < nx; i++) {
1403: PetscCall(VecGetSize(X[i], &Xng));
1404: PetscCall(VecGetLocalSize(X[i], &Xnl));
1405: PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1406: PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1407: shift += Xng;
1408: }
1409: /* create the concatenated vector */
1410: PetscCall(VecCreate(comm, &Ytmp));
1411: PetscCall(VecSetType(Ytmp, vec_type));
1412: PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1413: PetscCall(VecSetUp(Ytmp));
1414: /* copy data from X array to Y and return */
1415: for (i = 0; i < nx; i++) {
1416: PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1417: PetscCall(VecCopy(X[i], Xtmp));
1418: PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1419: }
1420: *Y = Ytmp;
1421: if (x_is) {
1422: *x_is = is_tmp;
1423: } else {
1424: for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1425: PetscCall(PetscFree(is_tmp));
1426: }
1427: }
1428: PetscFunctionReturn(PETSC_SUCCESS);
1429: }
1431: /* A helper function for VecGetSubVector to check if we can implement it with no-copy (i.e. the subvector shares
1432: memory with the original vector), and the block size of the subvector.
1434: Input Parameters:
1435: + X - the original vector
1436: - is - the index set of the subvector
1438: Output Parameters:
1439: + contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
1440: . start - start of contiguous block, as an offset from the start of the ownership range of the original vector
1441: - blocksize - the block size of the subvector
1443: */
1444: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1445: {
1446: PetscInt gstart, gend, lstart;
1447: PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1448: PetscInt n, N, ibs, vbs, bs = -1;
1450: PetscFunctionBegin;
1451: PetscCall(ISGetLocalSize(is, &n));
1452: PetscCall(ISGetSize(is, &N));
1453: PetscCall(ISGetBlockSize(is, &ibs));
1454: PetscCall(VecGetBlockSize(X, &vbs));
1455: PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1456: PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1457: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1458: if (ibs > 1) {
1459: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1460: bs = ibs;
1461: } else {
1462: if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1463: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1464: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1465: }
1467: *contig = red[0];
1468: *start = lstart;
1469: *blocksize = bs;
1470: PetscFunctionReturn(PETSC_SUCCESS);
1471: }
1473: /* A helper function for VecGetSubVector, to be used when we have to build a standalone subvector through VecScatter
1475: Input Parameters:
1476: + X - the original vector
1477: . is - the index set of the subvector
1478: - bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()
1480: Output Parameter:
1481: . Z - the subvector, which will compose the VecScatter context on output
1482: */
1483: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1484: {
1485: PetscInt n, N;
1486: VecScatter vscat;
1487: Vec Y;
1489: PetscFunctionBegin;
1490: PetscCall(ISGetLocalSize(is, &n));
1491: PetscCall(ISGetSize(is, &N));
1492: PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1493: PetscCall(VecSetSizes(Y, n, N));
1494: PetscCall(VecSetBlockSize(Y, bs));
1495: PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1496: PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1497: PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1498: PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1499: PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1500: PetscCall(VecScatterDestroy(&vscat));
1501: *Z = Y;
1502: PetscFunctionReturn(PETSC_SUCCESS);
1503: }
1505: /*@
1506: VecGetSubVector - Gets a vector representing part of another vector
1508: Collective
1510: Input Parameters:
1511: + X - vector from which to extract a subvector
1512: - is - index set representing portion of `X` to extract
1514: Output Parameter:
1515: . Y - subvector corresponding to `is`
1517: Level: advanced
1519: Notes:
1520: The subvector `Y` should be returned with `VecRestoreSubVector()`.
1521: `X` and must be defined on the same communicator
1523: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1524: modifying the subvector. Other non-overlapping subvectors can still be obtained from X using this function.
1526: 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`.
1528: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1529: @*/
1530: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1531: {
1532: Vec Z;
1534: PetscFunctionBegin;
1537: PetscCheckSameComm(X, 1, is, 2);
1538: PetscAssertPointer(Y, 3);
1539: if (X->ops->getsubvector) {
1540: PetscUseTypeMethod(X, getsubvector, is, &Z);
1541: } else { /* Default implementation currently does no caching */
1542: PetscBool contig;
1543: PetscInt n, N, start, bs;
1545: PetscCall(ISGetLocalSize(is, &n));
1546: PetscCall(ISGetSize(is, &N));
1547: PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1548: if (contig) { /* We can do a no-copy implementation */
1549: const PetscScalar *x;
1550: PetscInt state = 0;
1551: PetscBool isstd, iscuda, iship;
1553: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1554: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1555: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1556: if (iscuda) {
1557: #if defined(PETSC_HAVE_CUDA)
1558: const PetscScalar *x_d;
1559: PetscMPIInt size;
1560: PetscOffloadMask flg;
1562: PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1563: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1564: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1565: if (x) x += start;
1566: if (x_d) x_d += start;
1567: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1568: if (size == 1) {
1569: PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1570: } else {
1571: PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1572: }
1573: Z->offloadmask = flg;
1574: #endif
1575: } else if (iship) {
1576: #if defined(PETSC_HAVE_HIP)
1577: const PetscScalar *x_d;
1578: PetscMPIInt size;
1579: PetscOffloadMask flg;
1581: PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1582: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1583: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1584: if (x) x += start;
1585: if (x_d) x_d += start;
1586: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1587: if (size == 1) {
1588: PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1589: } else {
1590: PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1591: }
1592: Z->offloadmask = flg;
1593: #endif
1594: } else if (isstd) {
1595: PetscMPIInt size;
1597: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1598: PetscCall(VecGetArrayRead(X, &x));
1599: if (x) x += start;
1600: if (size == 1) {
1601: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1602: } else {
1603: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1604: }
1605: PetscCall(VecRestoreArrayRead(X, &x));
1606: } else { /* default implementation: use place array */
1607: PetscCall(VecGetArrayRead(X, &x));
1608: PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1609: PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1610: PetscCall(VecSetSizes(Z, n, N));
1611: PetscCall(VecSetBlockSize(Z, bs));
1612: PetscCall(VecPlaceArray(Z, x ? x + start : NULL));
1613: PetscCall(VecRestoreArrayRead(X, &x));
1614: }
1616: /* this is relevant only in debug mode */
1617: PetscCall(VecLockGet(X, &state));
1618: if (state) PetscCall(VecLockReadPush(Z));
1619: Z->ops->placearray = NULL;
1620: Z->ops->replacearray = NULL;
1621: } else { /* Have to create a scatter and do a copy */
1622: PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1623: }
1624: }
1625: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1626: if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1627: PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1628: *Y = Z;
1629: PetscFunctionReturn(PETSC_SUCCESS);
1630: }
1632: /*@
1633: VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`
1635: Collective
1637: Input Parameters:
1638: + X - vector from which subvector was obtained
1639: . is - index set representing the subset of `X`
1640: - Y - subvector being restored
1642: Level: advanced
1644: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1645: @*/
1646: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1647: {
1648: PETSC_UNUSED PetscObjectState dummystate = 0;
1649: PetscBool unchanged;
1651: PetscFunctionBegin;
1654: PetscCheckSameComm(X, 1, is, 2);
1655: PetscAssertPointer(Y, 3);
1658: if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1659: else {
1660: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1661: if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1662: VecScatter scatter;
1663: PetscInt state;
1665: PetscCall(VecLockGet(X, &state));
1666: PetscCheck(state == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec X is locked for read-only or read/write access");
1668: PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1669: if (scatter) {
1670: PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1671: PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1672: } else {
1673: PetscBool iscuda, iship;
1674: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1675: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1677: if (iscuda) {
1678: #if defined(PETSC_HAVE_CUDA)
1679: PetscOffloadMask ymask = (*Y)->offloadmask;
1681: /* The offloadmask of X dictates where to move memory
1682: If X GPU data is valid, then move Y data on GPU if needed
1683: Otherwise, move back to the CPU */
1684: switch (X->offloadmask) {
1685: case PETSC_OFFLOAD_BOTH:
1686: if (ymask == PETSC_OFFLOAD_CPU) {
1687: PetscCall(VecCUDAResetArray(*Y));
1688: } else if (ymask == PETSC_OFFLOAD_GPU) {
1689: X->offloadmask = PETSC_OFFLOAD_GPU;
1690: }
1691: break;
1692: case PETSC_OFFLOAD_GPU:
1693: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1694: break;
1695: case PETSC_OFFLOAD_CPU:
1696: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1697: break;
1698: case PETSC_OFFLOAD_UNALLOCATED:
1699: case PETSC_OFFLOAD_KOKKOS:
1700: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1701: }
1702: #endif
1703: } else if (iship) {
1704: #if defined(PETSC_HAVE_HIP)
1705: PetscOffloadMask ymask = (*Y)->offloadmask;
1707: /* The offloadmask of X dictates where to move memory
1708: If X GPU data is valid, then move Y data on GPU if needed
1709: Otherwise, move back to the CPU */
1710: switch (X->offloadmask) {
1711: case PETSC_OFFLOAD_BOTH:
1712: if (ymask == PETSC_OFFLOAD_CPU) {
1713: PetscCall(VecHIPResetArray(*Y));
1714: } else if (ymask == PETSC_OFFLOAD_GPU) {
1715: X->offloadmask = PETSC_OFFLOAD_GPU;
1716: }
1717: break;
1718: case PETSC_OFFLOAD_GPU:
1719: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1720: break;
1721: case PETSC_OFFLOAD_CPU:
1722: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1723: break;
1724: case PETSC_OFFLOAD_UNALLOCATED:
1725: case PETSC_OFFLOAD_KOKKOS:
1726: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1727: }
1728: #endif
1729: } else {
1730: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1731: PetscCall(VecResetArray(*Y));
1732: }
1733: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1734: }
1735: }
1736: }
1737: PetscCall(VecDestroy(Y));
1738: PetscFunctionReturn(PETSC_SUCCESS);
1739: }
1741: /*@
1742: VecCreateLocalVector - Creates a vector object suitable for use with `VecGetLocalVector()` and friends. You must call `VecDestroy()` when the
1743: vector is no longer needed.
1745: Not Collective.
1747: Input Parameter:
1748: . v - The vector for which the local vector is desired.
1750: Output Parameter:
1751: . w - Upon exit this contains the local vector.
1753: Level: beginner
1755: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1756: @*/
1757: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1758: {
1759: PetscMPIInt size;
1761: PetscFunctionBegin;
1763: PetscAssertPointer(w, 2);
1764: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1765: if (size == 1) PetscCall(VecDuplicate(v, w));
1766: else if (v->ops->createlocalvector) PetscUseTypeMethod(v, createlocalvector, w);
1767: else {
1768: VecType type;
1769: PetscInt n;
1771: PetscCall(VecCreate(PETSC_COMM_SELF, w));
1772: PetscCall(VecGetLocalSize(v, &n));
1773: PetscCall(VecSetSizes(*w, n, n));
1774: PetscCall(VecGetBlockSize(v, &n));
1775: PetscCall(VecSetBlockSize(*w, n));
1776: PetscCall(VecGetType(v, &type));
1777: PetscCall(VecSetType(*w, type));
1778: }
1779: PetscFunctionReturn(PETSC_SUCCESS);
1780: }
1782: /*@
1783: VecGetLocalVectorRead - Maps the local portion of a vector into a
1784: vector.
1786: Not Collective.
1788: Input Parameter:
1789: . v - The vector for which the local vector is desired.
1791: Output Parameter:
1792: . w - Upon exit this contains the local vector.
1794: Level: beginner
1796: Notes:
1797: You must call `VecRestoreLocalVectorRead()` when the local
1798: vector is no longer needed.
1800: This function is similar to `VecGetArrayRead()` which maps the local
1801: portion into a raw pointer. `VecGetLocalVectorRead()` is usually
1802: almost as efficient as `VecGetArrayRead()` but in certain circumstances
1803: `VecGetLocalVectorRead()` can be much more efficient than
1804: `VecGetArrayRead()`. This is because the construction of a contiguous
1805: array representing the vector data required by `VecGetArrayRead()` can
1806: be an expensive operation for certain vector types. For example, for
1807: GPU vectors `VecGetArrayRead()` requires that the data between device
1808: and host is synchronized.
1810: Unlike `VecGetLocalVector()`, this routine is not collective and
1811: preserves cached information.
1813: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1814: @*/
1815: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1816: {
1817: PetscFunctionBegin;
1820: VecCheckSameLocalSize(v, 1, w, 2);
1821: if (v->ops->getlocalvectorread) {
1822: PetscUseTypeMethod(v, getlocalvectorread, w);
1823: } else {
1824: PetscScalar *a;
1826: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1827: PetscCall(VecPlaceArray(w, a));
1828: }
1829: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1830: PetscCall(VecLockReadPush(v));
1831: PetscCall(VecLockReadPush(w));
1832: PetscFunctionReturn(PETSC_SUCCESS);
1833: }
1835: /*@
1836: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1837: previously mapped into a vector using `VecGetLocalVectorRead()`.
1839: Not Collective.
1841: Input Parameters:
1842: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
1843: - w - The vector into which the local portion of `v` was mapped.
1845: Level: beginner
1847: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1848: @*/
1849: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1850: {
1851: PetscFunctionBegin;
1854: if (v->ops->restorelocalvectorread) {
1855: PetscUseTypeMethod(v, restorelocalvectorread, w);
1856: } else {
1857: const PetscScalar *a;
1859: PetscCall(VecGetArrayRead(w, &a));
1860: PetscCall(VecRestoreArrayRead(v, &a));
1861: PetscCall(VecResetArray(w));
1862: }
1863: PetscCall(VecLockReadPop(v));
1864: PetscCall(VecLockReadPop(w));
1865: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1866: PetscFunctionReturn(PETSC_SUCCESS);
1867: }
1869: /*@
1870: VecGetLocalVector - Maps the local portion of a vector into a
1871: vector.
1873: Collective
1875: Input Parameter:
1876: . v - The vector for which the local vector is desired.
1878: Output Parameter:
1879: . w - Upon exit this contains the local vector.
1881: Level: beginner
1883: Notes:
1884: You must call `VecRestoreLocalVector()` when the local
1885: vector is no longer needed.
1887: This function is similar to `VecGetArray()` which maps the local
1888: portion into a raw pointer. `VecGetLocalVector()` is usually about as
1889: efficient as `VecGetArray()` but in certain circumstances
1890: `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1891: This is because the construction of a contiguous array representing
1892: the vector data required by `VecGetArray()` can be an expensive
1893: operation for certain vector types. For example, for GPU vectors
1894: `VecGetArray()` requires that the data between device and host is
1895: synchronized.
1897: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1898: @*/
1899: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1900: {
1901: PetscFunctionBegin;
1904: VecCheckSameLocalSize(v, 1, w, 2);
1905: if (v->ops->getlocalvector) {
1906: PetscUseTypeMethod(v, getlocalvector, w);
1907: } else {
1908: PetscScalar *a;
1910: PetscCall(VecGetArray(v, &a));
1911: PetscCall(VecPlaceArray(w, a));
1912: }
1913: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1914: PetscFunctionReturn(PETSC_SUCCESS);
1915: }
1917: /*@
1918: VecRestoreLocalVector - Unmaps the local portion of a vector
1919: previously mapped into a vector using `VecGetLocalVector()`.
1921: Logically Collective.
1923: Input Parameters:
1924: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
1925: - w - The vector into which the local portion of `v` was mapped.
1927: Level: beginner
1929: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1930: @*/
1931: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1932: {
1933: PetscFunctionBegin;
1936: if (v->ops->restorelocalvector) {
1937: PetscUseTypeMethod(v, restorelocalvector, w);
1938: } else {
1939: PetscScalar *a;
1940: PetscCall(VecGetArray(w, &a));
1941: PetscCall(VecRestoreArray(v, &a));
1942: PetscCall(VecResetArray(w));
1943: }
1944: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1945: PetscCall(PetscObjectStateIncrease((PetscObject)v));
1946: PetscFunctionReturn(PETSC_SUCCESS);
1947: }
1949: /*@C
1950: VecGetArray - Returns a pointer to a contiguous array that contains this
1951: MPI processes's portion of the vector data
1953: Logically Collective
1955: Input Parameter:
1956: . x - the vector
1958: Output Parameter:
1959: . a - location to put pointer to the array
1961: Level: beginner
1963: Notes:
1964: For the standard PETSc vectors, `VecGetArray()` returns a pointer to the local data array and
1965: does not use any copies. If the underlying vector data is not stored in a contiguous array
1966: this routine will copy the data to a contiguous array and return a pointer to that. You MUST
1967: call `VecRestoreArray()` when you no longer need access to the array.
1969: Fortran Notes:
1970: `VecGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayF90()`
1972: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
1973: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
1974: @*/
1975: PetscErrorCode VecGetArray(Vec x, PetscScalar **a)
1976: {
1977: PetscFunctionBegin;
1979: PetscCall(VecSetErrorIfLocked(x, 1));
1980: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
1981: PetscUseTypeMethod(x, getarray, a);
1982: } else if (x->petscnative) { /* VECSTANDARD */
1983: *a = *((PetscScalar **)x->data);
1984: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
1985: PetscFunctionReturn(PETSC_SUCCESS);
1986: }
1988: /*@C
1989: VecRestoreArray - Restores a vector after `VecGetArray()` has been called and the array is no longer needed
1991: Logically Collective
1993: Input Parameters:
1994: + x - the vector
1995: - a - location of pointer to array obtained from `VecGetArray()`
1997: Level: beginner
1999: Fortran Notes:
2000: `VecRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayF90()`
2002: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2003: `VecGetArrayPair()`, `VecRestoreArrayPair()`
2004: @*/
2005: PetscErrorCode VecRestoreArray(Vec x, PetscScalar **a)
2006: {
2007: PetscFunctionBegin;
2009: if (a) PetscAssertPointer(a, 2);
2010: if (x->ops->restorearray) {
2011: PetscUseTypeMethod(x, restorearray, a);
2012: } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
2013: if (a) *a = NULL;
2014: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2015: PetscFunctionReturn(PETSC_SUCCESS);
2016: }
2017: /*@C
2018: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
2020: Not Collective
2022: Input Parameter:
2023: . x - the vector
2025: Output Parameter:
2026: . a - the array
2028: Level: beginner
2030: Notes:
2031: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2033: Unlike `VecGetArray()`, preserves cached information like vector norms.
2035: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
2036: implementations may require a copy, but such implementations should cache the contiguous representation so that
2037: only one copy is performed when this routine is called multiple times in sequence.
2039: Fortran Notes:
2040: `VecGetArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayReadF90()`
2042: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2043: @*/
2044: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar **a)
2045: {
2046: PetscFunctionBegin;
2048: PetscAssertPointer(a, 2);
2049: if (x->ops->getarrayread) {
2050: PetscUseTypeMethod(x, getarrayread, a);
2051: } else if (x->ops->getarray) {
2052: PetscObjectState state;
2054: /* VECNEST, VECCUDA, VECKOKKOS etc */
2055: // x->ops->getarray may bump the object state, but since we know this is a read-only get
2056: // we can just undo that
2057: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2058: PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
2059: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2060: } else if (x->petscnative) {
2061: /* VECSTANDARD */
2062: *a = *((PetscScalar **)x->data);
2063: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
2064: PetscFunctionReturn(PETSC_SUCCESS);
2065: }
2067: /*@C
2068: VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`
2070: Not Collective
2072: Input Parameters:
2073: + x - the vector
2074: - a - the array
2076: Level: beginner
2078: Fortran Notes:
2079: `VecRestoreArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayReadF90()`
2081: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2082: @*/
2083: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar **a)
2084: {
2085: PetscFunctionBegin;
2087: if (a) PetscAssertPointer(a, 2);
2088: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
2089: /* nothing */
2090: } else if (x->ops->restorearrayread) { /* VECNEST */
2091: PetscUseTypeMethod(x, restorearrayread, a);
2092: } else { /* No one? */
2093: PetscObjectState state;
2095: // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2096: // we can just undo that
2097: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2098: PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2099: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2100: }
2101: if (a) *a = NULL;
2102: PetscFunctionReturn(PETSC_SUCCESS);
2103: }
2105: /*@C
2106: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
2107: MPI processes's portion of the vector data.
2109: Logically Collective
2111: Input Parameter:
2112: . x - the vector
2114: Output Parameter:
2115: . a - location to put pointer to the array
2117: Level: intermediate
2119: Note:
2120: The values in this array are NOT valid, the caller of this routine is responsible for putting
2121: values into the array; any values it does not set will be invalid.
2123: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2125: For vectors associated with GPUs, the host and device vectors are not synchronized before
2126: giving access. If you need correct values in the array use `VecGetArray()`
2128: Fortran Notes:
2129: `VecGetArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayWriteF90()`
2131: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteF90()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
2132: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
2133: @*/
2134: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar **a)
2135: {
2136: PetscFunctionBegin;
2138: PetscAssertPointer(a, 2);
2139: PetscCall(VecSetErrorIfLocked(x, 1));
2140: if (x->ops->getarraywrite) {
2141: PetscUseTypeMethod(x, getarraywrite, a);
2142: } else {
2143: PetscCall(VecGetArray(x, a));
2144: }
2145: PetscFunctionReturn(PETSC_SUCCESS);
2146: }
2148: /*@C
2149: VecRestoreArrayWrite - Restores a vector after `VecGetArrayWrite()` has been called.
2151: Logically Collective
2153: Input Parameters:
2154: + x - the vector
2155: - a - location of pointer to array obtained from `VecGetArray()`
2157: Level: beginner
2159: Fortran Notes:
2160: `VecRestoreArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayWriteF90()`
2162: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteF90()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2163: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2164: @*/
2165: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar **a)
2166: {
2167: PetscFunctionBegin;
2169: if (a) PetscAssertPointer(a, 2);
2170: if (x->ops->restorearraywrite) {
2171: PetscUseTypeMethod(x, restorearraywrite, a);
2172: } else if (x->ops->restorearray) {
2173: PetscUseTypeMethod(x, restorearray, a);
2174: }
2175: if (a) *a = NULL;
2176: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2177: PetscFunctionReturn(PETSC_SUCCESS);
2178: }
2180: /*@C
2181: VecGetArrays - Returns a pointer to the arrays in a set of vectors
2182: that were created by a call to `VecDuplicateVecs()`.
2184: Logically Collective; No Fortran Support
2186: Input Parameters:
2187: + x - the vectors
2188: - n - the number of vectors
2190: Output Parameter:
2191: . a - location to put pointer to the array
2193: Level: intermediate
2195: Note:
2196: You MUST call `VecRestoreArrays()` when you no longer need access to the arrays.
2198: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2199: @*/
2200: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2201: {
2202: PetscInt i;
2203: PetscScalar **q;
2205: PetscFunctionBegin;
2206: PetscAssertPointer(x, 1);
2208: PetscAssertPointer(a, 3);
2209: PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2210: PetscCall(PetscMalloc1(n, &q));
2211: for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2212: *a = q;
2213: PetscFunctionReturn(PETSC_SUCCESS);
2214: }
2216: /*@C
2217: VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2218: has been called.
2220: Logically Collective; No Fortran Support
2222: Input Parameters:
2223: + x - the vector
2224: . n - the number of vectors
2225: - a - location of pointer to arrays obtained from `VecGetArrays()`
2227: Notes:
2228: For regular PETSc vectors this routine does not involve any copies. For
2229: any special vectors that do not store local vector data in a contiguous
2230: array, this routine will copy the data back into the underlying
2231: vector data structure from the arrays obtained with `VecGetArrays()`.
2233: Level: intermediate
2235: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2236: @*/
2237: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2238: {
2239: PetscInt i;
2240: PetscScalar **q = *a;
2242: PetscFunctionBegin;
2243: PetscAssertPointer(x, 1);
2245: PetscAssertPointer(a, 3);
2247: for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2248: PetscCall(PetscFree(q));
2249: PetscFunctionReturn(PETSC_SUCCESS);
2250: }
2252: /*@C
2253: VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g.,
2254: `VECCUDA`), the returned pointer will be a device pointer to the device memory that contains
2255: this MPI processes's portion of the vector data.
2257: Logically Collective; No Fortran Support
2259: Input Parameter:
2260: . x - the vector
2262: Output Parameters:
2263: + a - location to put pointer to the array
2264: - mtype - memory type of the array
2266: Level: beginner
2268: Note:
2269: Device data is guaranteed to have the latest value. Otherwise, when this is a host vector
2270: (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host
2271: pointer.
2273: For `VECKOKKOS`, if Kokkos is configured without device (e.g., use serial or openmp), per
2274: this function, the vector works like `VECSEQ`/`VECMPI`; otherwise, it works like `VECCUDA` or
2275: `VECHIP` etc.
2277: Use `VecRestoreArrayAndMemType()` when the array access is no longer needed.
2279: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`,
2280: `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2281: @*/
2282: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2283: {
2284: PetscFunctionBegin;
2287: PetscAssertPointer(a, 2);
2288: if (mtype) PetscAssertPointer(mtype, 3);
2289: PetscCall(VecSetErrorIfLocked(x, 1));
2290: if (x->ops->getarrayandmemtype) {
2291: /* VECCUDA, VECKOKKOS etc */
2292: PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2293: } else {
2294: /* VECSTANDARD, VECNEST, VECVIENNACL */
2295: PetscCall(VecGetArray(x, a));
2296: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2297: }
2298: PetscFunctionReturn(PETSC_SUCCESS);
2299: }
2301: /*@C
2302: VecRestoreArrayAndMemType - Restores a vector after `VecGetArrayAndMemType()` has been called.
2304: Logically Collective; No Fortran Support
2306: Input Parameters:
2307: + x - the vector
2308: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`
2310: Level: beginner
2312: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`,
2313: `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2314: @*/
2315: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar **a)
2316: {
2317: PetscFunctionBegin;
2320: if (a) PetscAssertPointer(a, 2);
2321: if (x->ops->restorearrayandmemtype) {
2322: /* VECCUDA, VECKOKKOS etc */
2323: PetscUseTypeMethod(x, restorearrayandmemtype, a);
2324: } else {
2325: /* VECNEST, VECVIENNACL */
2326: PetscCall(VecRestoreArray(x, a));
2327: } /* VECSTANDARD does nothing */
2328: if (a) *a = NULL;
2329: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2330: PetscFunctionReturn(PETSC_SUCCESS);
2331: }
2333: /*@C
2334: VecGetArrayReadAndMemType - Like `VecGetArrayRead()`, but if the input vector is a device vector, it will return a read-only device pointer.
2335: The returned pointer is guaranteed to point to up-to-date data. For host vectors, it functions as `VecGetArrayRead()`.
2337: Not Collective; No Fortran Support
2339: Input Parameter:
2340: . x - the vector
2342: Output Parameters:
2343: + a - the array
2344: - mtype - memory type of the array
2346: Level: beginner
2348: Notes:
2349: The array must be returned using a matching call to `VecRestoreArrayReadAndMemType()`.
2351: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2352: @*/
2353: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar **a, PetscMemType *mtype)
2354: {
2355: PetscFunctionBegin;
2358: PetscAssertPointer(a, 2);
2359: if (mtype) PetscAssertPointer(mtype, 3);
2360: if (x->ops->getarrayreadandmemtype) {
2361: /* VECCUDA/VECHIP though they are also petscnative */
2362: PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2363: } else if (x->ops->getarrayandmemtype) {
2364: /* VECKOKKOS */
2365: PetscObjectState state;
2367: // see VecGetArrayRead() for why
2368: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2369: PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2370: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2371: } else {
2372: PetscCall(VecGetArrayRead(x, a));
2373: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2374: }
2375: PetscFunctionReturn(PETSC_SUCCESS);
2376: }
2378: /*@C
2379: VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`
2381: Not Collective; No Fortran Support
2383: Input Parameters:
2384: + x - the vector
2385: - a - the array
2387: Level: beginner
2389: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2390: @*/
2391: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar **a)
2392: {
2393: PetscFunctionBegin;
2396: if (a) PetscAssertPointer(a, 2);
2397: if (x->ops->restorearrayreadandmemtype) {
2398: /* VECCUDA/VECHIP */
2399: PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2400: } else if (!x->petscnative) {
2401: /* VECNEST */
2402: PetscCall(VecRestoreArrayRead(x, a));
2403: }
2404: if (a) *a = NULL;
2405: PetscFunctionReturn(PETSC_SUCCESS);
2406: }
2408: /*@C
2409: VecGetArrayWriteAndMemType - Like `VecGetArrayWrite()`, but if this is a device vector it will always return
2410: a device pointer to the device memory that contains this processor's portion of the vector data.
2412: Logically Collective; No Fortran Support
2414: Input Parameter:
2415: . x - the vector
2417: Output Parameters:
2418: + a - the array
2419: - mtype - memory type of the array
2421: Level: beginner
2423: Note:
2424: The array must be returned using a matching call to `VecRestoreArrayWriteAndMemType()`, where it will label the device memory as most recent.
2426: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2427: @*/
2428: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2429: {
2430: PetscFunctionBegin;
2433: PetscCall(VecSetErrorIfLocked(x, 1));
2434: PetscAssertPointer(a, 2);
2435: if (mtype) PetscAssertPointer(mtype, 3);
2436: if (x->ops->getarraywriteandmemtype) {
2437: /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2438: PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2439: } else if (x->ops->getarrayandmemtype) {
2440: PetscCall(VecGetArrayAndMemType(x, a, mtype));
2441: } else {
2442: /* VECNEST, VECVIENNACL */
2443: PetscCall(VecGetArrayWrite(x, a));
2444: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2445: }
2446: PetscFunctionReturn(PETSC_SUCCESS);
2447: }
2449: /*@C
2450: VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`
2452: Logically Collective; No Fortran Support
2454: Input Parameters:
2455: + x - the vector
2456: - a - the array
2458: Level: beginner
2460: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2461: @*/
2462: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar **a)
2463: {
2464: PetscFunctionBegin;
2467: PetscCall(VecSetErrorIfLocked(x, 1));
2468: if (a) PetscAssertPointer(a, 2);
2469: if (x->ops->restorearraywriteandmemtype) {
2470: /* VECCUDA/VECHIP */
2471: PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2472: PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2473: } else if (x->ops->restorearrayandmemtype) {
2474: PetscCall(VecRestoreArrayAndMemType(x, a));
2475: } else {
2476: PetscCall(VecRestoreArray(x, a));
2477: }
2478: if (a) *a = NULL;
2479: PetscFunctionReturn(PETSC_SUCCESS);
2480: }
2482: /*@
2483: VecPlaceArray - Allows one to replace the array in a vector with an
2484: array provided by the user. This is useful to avoid copying an array
2485: into a vector.
2487: Logically Collective; No Fortran Support
2489: Input Parameters:
2490: + vec - the vector
2491: - array - the array
2493: Level: developer
2495: Notes:
2496: Use `VecReplaceArray()` instead to permanently replace the array
2498: You can return to the original array with a call to `VecResetArray()`. `vec` does not take
2499: ownership of `array` in any way.
2501: The user must free `array` themselves but be careful not to
2502: do so before the vector has either been destroyed, had its original array restored with
2503: `VecResetArray()` or permanently replaced with `VecReplaceArray()`.
2505: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2506: @*/
2507: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2508: {
2509: PetscFunctionBegin;
2512: if (array) PetscAssertPointer(array, 2);
2513: PetscUseTypeMethod(vec, placearray, array);
2514: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2515: PetscFunctionReturn(PETSC_SUCCESS);
2516: }
2518: /*@C
2519: VecReplaceArray - Allows one to replace the array in a vector with an
2520: array provided by the user. This is useful to avoid copying an array
2521: into a vector.
2523: Logically Collective; No Fortran Support
2525: Input Parameters:
2526: + vec - the vector
2527: - array - the array
2529: Level: developer
2531: Notes:
2532: This permanently replaces the array and frees the memory associated
2533: with the old array. Use `VecPlaceArray()` to temporarily replace the array.
2535: The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
2536: freed by the user. It will be freed when the vector is destroyed.
2538: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2539: @*/
2540: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2541: {
2542: PetscFunctionBegin;
2545: PetscUseTypeMethod(vec, replacearray, array);
2546: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2547: PetscFunctionReturn(PETSC_SUCCESS);
2548: }
2550: /*MC
2551: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
2552: and makes them accessible via a Fortran pointer.
2554: Synopsis:
2555: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
2557: Collective
2559: Input Parameters:
2560: + x - a vector to mimic
2561: - n - the number of vectors to obtain
2563: Output Parameters:
2564: + y - Fortran pointer to the array of vectors
2565: - ierr - error code
2567: Example of Usage:
2568: .vb
2569: #include <petsc/finclude/petscvec.h>
2570: use petscvec
2572: Vec x
2573: Vec, pointer :: y(:)
2574: ....
2575: call VecDuplicateVecsF90(x,2,y,ierr)
2576: call VecSet(y(2),alpha,ierr)
2577: call VecSet(y(2),alpha,ierr)
2578: ....
2579: call VecDestroyVecsF90(2,y,ierr)
2580: .ve
2582: Level: beginner
2584: Note:
2585: Use `VecDestroyVecsF90()` to free the space.
2587: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecsF90()`, `VecDuplicateVecs()`
2588: M*/
2590: /*MC
2591: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2592: `VecGetArrayF90()`.
2594: Synopsis:
2595: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2597: Logically Collective
2599: Input Parameters:
2600: + x - vector
2601: - xx_v - the Fortran pointer to the array
2603: Output Parameter:
2604: . ierr - error code
2606: Example of Usage:
2607: .vb
2608: #include <petsc/finclude/petscvec.h>
2609: use petscvec
2611: PetscScalar, pointer :: xx_v(:)
2612: ....
2613: call VecGetArrayF90(x,xx_v,ierr)
2614: xx_v(3) = a
2615: call VecRestoreArrayF90(x,xx_v,ierr)
2616: .ve
2618: Level: beginner
2620: .seealso: [](ch_vectors), `Vec`, `VecGetArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrayReadF90()`
2621: M*/
2623: /*MC
2624: VecDestroyVecsF90 - Frees a block of vectors obtained with `VecDuplicateVecsF90()`.
2626: Synopsis:
2627: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
2629: Collective
2631: Input Parameters:
2632: + n - the number of vectors previously obtained
2633: - x - pointer to array of vector pointers
2635: Output Parameter:
2636: . ierr - error code
2638: Level: beginner
2640: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecs()`, `VecDuplicateVecsF90()`
2641: M*/
2643: /*MC
2644: VecGetArrayF90 - Accesses a vector array from Fortran. For default PETSc
2645: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2646: this routine is implementation dependent. You MUST call `VecRestoreArrayF90()`
2647: when you no longer need access to the array.
2649: Synopsis:
2650: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2652: Logically Collective
2654: Input Parameter:
2655: . x - vector
2657: Output Parameters:
2658: + xx_v - the Fortran pointer to the array
2659: - ierr - error code
2661: Example of Usage:
2662: .vb
2663: #include <petsc/finclude/petscvec.h>
2664: use petscvec
2666: PetscScalar, pointer :: xx_v(:)
2667: ....
2668: call VecGetArrayF90(x,xx_v,ierr)
2669: xx_v(3) = a
2670: call VecRestoreArrayF90(x,xx_v,ierr)
2671: .ve
2673: Level: beginner
2675: Note:
2676: If you ONLY intend to read entries from the array and not change any entries you should use `VecGetArrayReadF90()`.
2678: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayReadF90()`
2679: M*/
2681: /*MC
2682: VecGetArrayReadF90 - Accesses a read only array from Fortran. For default PETSc
2683: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2684: this routine is implementation dependent. You MUST call `VecRestoreArrayReadF90()`
2685: when you no longer need access to the array.
2687: Synopsis:
2688: VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2690: Logically Collective
2692: Input Parameter:
2693: . x - vector
2695: Output Parameters:
2696: + xx_v - the Fortran pointer to the array
2697: - ierr - error code
2699: Example of Usage:
2700: .vb
2701: #include <petsc/finclude/petscvec.h>
2702: use petscvec
2704: PetscScalar, pointer :: xx_v(:)
2705: ....
2706: call VecGetArrayReadF90(x,xx_v,ierr)
2707: a = xx_v(3)
2708: call VecRestoreArrayReadF90(x,xx_v,ierr)
2709: .ve
2711: Level: beginner
2713: Note:
2714: If you intend to write entries into the array you must use `VecGetArrayF90()`.
2716: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecGetArrayF90()`
2717: M*/
2719: /*MC
2720: VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2721: `VecGetArrayReadF90()`.
2723: Synopsis:
2724: VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2726: Logically Collective
2728: Input Parameters:
2729: + x - vector
2730: - xx_v - the Fortran pointer to the array
2732: Output Parameter:
2733: . ierr - error code
2735: Example of Usage:
2736: .vb
2737: #include <petsc/finclude/petscvec.h>
2738: use petscvec
2740: PetscScalar, pointer :: xx_v(:)
2741: ....
2742: call VecGetArrayReadF90(x,xx_v,ierr)
2743: a = xx_v(3)
2744: call VecRestoreArrayReadF90(x,xx_v,ierr)
2745: .ve
2747: Level: beginner
2749: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecRestoreArrayF90()`
2750: M*/
2752: /*@C
2753: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2754: processor's portion of the vector data. You MUST call `VecRestoreArray2d()`
2755: when you no longer need access to the array.
2757: Logically Collective
2759: Input Parameters:
2760: + x - the vector
2761: . m - first dimension of two dimensional array
2762: . n - second dimension of two dimensional array
2763: . mstart - first index you will use in first coordinate direction (often 0)
2764: - nstart - first index in the second coordinate direction (often 0)
2766: Output Parameter:
2767: . a - location to put pointer to the array
2769: Level: developer
2771: Notes:
2772: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2773: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2774: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2775: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2777: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2779: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2780: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2781: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2782: @*/
2783: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2784: {
2785: PetscInt i, N;
2786: PetscScalar *aa;
2788: PetscFunctionBegin;
2790: PetscAssertPointer(a, 6);
2792: PetscCall(VecGetLocalSize(x, &N));
2793: 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);
2794: PetscCall(VecGetArray(x, &aa));
2796: PetscCall(PetscMalloc1(m, a));
2797: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2798: *a -= mstart;
2799: PetscFunctionReturn(PETSC_SUCCESS);
2800: }
2802: /*@C
2803: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2804: processor's portion of the vector data. You MUST call `VecRestoreArray2dWrite()`
2805: when you no longer need access to the array.
2807: Logically Collective
2809: Input Parameters:
2810: + x - the vector
2811: . m - first dimension of two dimensional array
2812: . n - second dimension of two dimensional array
2813: . mstart - first index you will use in first coordinate direction (often 0)
2814: - nstart - first index in the second coordinate direction (often 0)
2816: Output Parameter:
2817: . a - location to put pointer to the array
2819: Level: developer
2821: Notes:
2822: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2823: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2824: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2825: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2827: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2829: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2830: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2831: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2832: @*/
2833: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2834: {
2835: PetscInt i, N;
2836: PetscScalar *aa;
2838: PetscFunctionBegin;
2840: PetscAssertPointer(a, 6);
2842: PetscCall(VecGetLocalSize(x, &N));
2843: 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);
2844: PetscCall(VecGetArrayWrite(x, &aa));
2846: PetscCall(PetscMalloc1(m, a));
2847: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2848: *a -= mstart;
2849: PetscFunctionReturn(PETSC_SUCCESS);
2850: }
2852: /*@C
2853: VecRestoreArray2d - Restores a vector after `VecGetArray2d()` has been called.
2855: Logically Collective
2857: Input Parameters:
2858: + x - the vector
2859: . m - first dimension of two dimensional array
2860: . n - second dimension of the two dimensional array
2861: . mstart - first index you will use in first coordinate direction (often 0)
2862: . nstart - first index in the second coordinate direction (often 0)
2863: - a - location of pointer to array obtained from `VecGetArray2d()`
2865: Level: developer
2867: Notes:
2868: For regular PETSc vectors this routine does not involve any copies. For
2869: any special vectors that do not store local vector data in a contiguous
2870: array, this routine will copy the data back into the underlying
2871: vector data structure from the array obtained with `VecGetArray()`.
2873: This routine actually zeros out the `a` pointer.
2875: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2876: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2877: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2878: @*/
2879: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2880: {
2881: void *dummy;
2883: PetscFunctionBegin;
2885: PetscAssertPointer(a, 6);
2887: dummy = (void *)(*a + mstart);
2888: PetscCall(PetscFree(dummy));
2889: PetscCall(VecRestoreArray(x, NULL));
2890: *a = NULL;
2891: PetscFunctionReturn(PETSC_SUCCESS);
2892: }
2894: /*@C
2895: VecRestoreArray2dWrite - Restores a vector after `VecGetArray2dWrite()` has been called.
2897: Logically Collective
2899: Input Parameters:
2900: + x - the vector
2901: . m - first dimension of two dimensional array
2902: . n - second dimension of the two dimensional array
2903: . mstart - first index you will use in first coordinate direction (often 0)
2904: . nstart - first index in the second coordinate direction (often 0)
2905: - a - location of pointer to array obtained from `VecGetArray2d()`
2907: Level: developer
2909: Notes:
2910: For regular PETSc vectors this routine does not involve any copies. For
2911: any special vectors that do not store local vector data in a contiguous
2912: array, this routine will copy the data back into the underlying
2913: vector data structure from the array obtained with `VecGetArray()`.
2915: This routine actually zeros out the `a` pointer.
2917: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2918: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2919: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2920: @*/
2921: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2922: {
2923: void *dummy;
2925: PetscFunctionBegin;
2927: PetscAssertPointer(a, 6);
2929: dummy = (void *)(*a + mstart);
2930: PetscCall(PetscFree(dummy));
2931: PetscCall(VecRestoreArrayWrite(x, NULL));
2932: PetscFunctionReturn(PETSC_SUCCESS);
2933: }
2935: /*@C
2936: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2937: processor's portion of the vector data. You MUST call `VecRestoreArray1d()`
2938: when you no longer need access to the array.
2940: Logically Collective
2942: Input Parameters:
2943: + x - the vector
2944: . m - first dimension of two dimensional array
2945: - mstart - first index you will use in first coordinate direction (often 0)
2947: Output Parameter:
2948: . a - location to put pointer to the array
2950: Level: developer
2952: Notes:
2953: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2954: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2955: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2957: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2959: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2960: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2961: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2962: @*/
2963: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2964: {
2965: PetscInt N;
2967: PetscFunctionBegin;
2969: PetscAssertPointer(a, 4);
2971: PetscCall(VecGetLocalSize(x, &N));
2972: 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);
2973: PetscCall(VecGetArray(x, a));
2974: *a -= mstart;
2975: PetscFunctionReturn(PETSC_SUCCESS);
2976: }
2978: /*@C
2979: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
2980: processor's portion of the vector data. You MUST call `VecRestoreArray1dWrite()`
2981: when you no longer need access to the array.
2983: Logically Collective
2985: Input Parameters:
2986: + x - the vector
2987: . m - first dimension of two dimensional array
2988: - mstart - first index you will use in first coordinate direction (often 0)
2990: Output Parameter:
2991: . a - location to put pointer to the array
2993: Level: developer
2995: Notes:
2996: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2997: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2998: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3000: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3002: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3003: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3004: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3005: @*/
3006: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3007: {
3008: PetscInt N;
3010: PetscFunctionBegin;
3012: PetscAssertPointer(a, 4);
3014: PetscCall(VecGetLocalSize(x, &N));
3015: 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);
3016: PetscCall(VecGetArrayWrite(x, a));
3017: *a -= mstart;
3018: PetscFunctionReturn(PETSC_SUCCESS);
3019: }
3021: /*@C
3022: VecRestoreArray1d - Restores a vector after `VecGetArray1d()` has been called.
3024: Logically Collective
3026: Input Parameters:
3027: + x - the vector
3028: . m - first dimension of two dimensional array
3029: . mstart - first index you will use in first coordinate direction (often 0)
3030: - a - location of pointer to array obtained from `VecGetArray1d()`
3032: Level: developer
3034: Notes:
3035: For regular PETSc vectors this routine does not involve any copies. For
3036: any special vectors that do not store local vector data in a contiguous
3037: array, this routine will copy the data back into the underlying
3038: vector data structure from the array obtained with `VecGetArray1d()`.
3040: This routine actually zeros out the `a` pointer.
3042: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3043: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3044: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3045: @*/
3046: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3047: {
3048: PetscFunctionBegin;
3051: PetscCall(VecRestoreArray(x, NULL));
3052: *a = NULL;
3053: PetscFunctionReturn(PETSC_SUCCESS);
3054: }
3056: /*@C
3057: VecRestoreArray1dWrite - Restores a vector after `VecGetArray1dWrite()` has been called.
3059: Logically Collective
3061: Input Parameters:
3062: + x - the vector
3063: . m - first dimension of two dimensional array
3064: . mstart - first index you will use in first coordinate direction (often 0)
3065: - a - location of pointer to array obtained from `VecGetArray1d()`
3067: Level: developer
3069: Notes:
3070: For regular PETSc vectors this routine does not involve any copies. For
3071: any special vectors that do not store local vector data in a contiguous
3072: array, this routine will copy the data back into the underlying
3073: vector data structure from the array obtained with `VecGetArray1d()`.
3075: This routine actually zeros out the `a` pointer.
3077: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3078: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3079: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3080: @*/
3081: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3082: {
3083: PetscFunctionBegin;
3086: PetscCall(VecRestoreArrayWrite(x, NULL));
3087: *a = NULL;
3088: PetscFunctionReturn(PETSC_SUCCESS);
3089: }
3091: /*@C
3092: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
3093: processor's portion of the vector data. You MUST call `VecRestoreArray3d()`
3094: when you no longer need access to the array.
3096: Logically Collective
3098: Input Parameters:
3099: + x - the vector
3100: . m - first dimension of three dimensional array
3101: . n - second dimension of three dimensional array
3102: . p - third dimension of three dimensional array
3103: . mstart - first index you will use in first coordinate direction (often 0)
3104: . nstart - first index in the second coordinate direction (often 0)
3105: - pstart - first index in the third coordinate direction (often 0)
3107: Output Parameter:
3108: . a - location to put pointer to the array
3110: Level: developer
3112: Notes:
3113: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3114: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3115: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3116: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3118: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3120: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3121: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
3122: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3123: @*/
3124: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3125: {
3126: PetscInt i, N, j;
3127: PetscScalar *aa, **b;
3129: PetscFunctionBegin;
3131: PetscAssertPointer(a, 8);
3133: PetscCall(VecGetLocalSize(x, &N));
3134: 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);
3135: PetscCall(VecGetArray(x, &aa));
3137: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3138: b = (PetscScalar **)((*a) + m);
3139: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3140: for (i = 0; i < m; i++)
3141: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3142: *a -= mstart;
3143: PetscFunctionReturn(PETSC_SUCCESS);
3144: }
3146: /*@C
3147: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
3148: processor's portion of the vector data. You MUST call `VecRestoreArray3dWrite()`
3149: when you no longer need access to the array.
3151: Logically Collective
3153: Input Parameters:
3154: + x - the vector
3155: . m - first dimension of three dimensional array
3156: . n - second dimension of three dimensional array
3157: . p - third dimension of three 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)
3162: Output Parameter:
3163: . a - location to put pointer to the array
3165: Level: developer
3167: Notes:
3168: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3169: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3170: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3171: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3173: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3175: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3176: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3177: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3178: @*/
3179: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3180: {
3181: PetscInt i, N, j;
3182: PetscScalar *aa, **b;
3184: PetscFunctionBegin;
3186: PetscAssertPointer(a, 8);
3188: PetscCall(VecGetLocalSize(x, &N));
3189: 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);
3190: PetscCall(VecGetArrayWrite(x, &aa));
3192: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3193: b = (PetscScalar **)((*a) + m);
3194: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3195: for (i = 0; i < m; i++)
3196: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3198: *a -= mstart;
3199: PetscFunctionReturn(PETSC_SUCCESS);
3200: }
3202: /*@C
3203: VecRestoreArray3d - Restores a vector after `VecGetArray3d()` has been called.
3205: Logically Collective
3207: Input Parameters:
3208: + x - the vector
3209: . m - first dimension of three dimensional array
3210: . n - second dimension of the three dimensional array
3211: . p - third dimension of the three dimensional array
3212: . mstart - first index you will use in first coordinate direction (often 0)
3213: . nstart - first index in the second coordinate direction (often 0)
3214: . pstart - first index in the third coordinate direction (often 0)
3215: - a - location of pointer to array obtained from VecGetArray3d()
3217: Level: developer
3219: Notes:
3220: For regular PETSc vectors this routine does not involve any copies. For
3221: any special vectors that do not store local vector data in a contiguous
3222: array, this routine will copy the data back into the underlying
3223: vector data structure from the array obtained with `VecGetArray()`.
3225: This routine actually zeros out the `a` pointer.
3227: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3228: `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3229: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3230: @*/
3231: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3232: {
3233: void *dummy;
3235: PetscFunctionBegin;
3237: PetscAssertPointer(a, 8);
3239: dummy = (void *)(*a + mstart);
3240: PetscCall(PetscFree(dummy));
3241: PetscCall(VecRestoreArray(x, NULL));
3242: *a = NULL;
3243: PetscFunctionReturn(PETSC_SUCCESS);
3244: }
3246: /*@C
3247: VecRestoreArray3dWrite - Restores a vector after `VecGetArray3dWrite()` has been called.
3249: Logically Collective
3251: Input Parameters:
3252: + x - the vector
3253: . m - first dimension of three dimensional array
3254: . n - second dimension of the three dimensional array
3255: . p - third dimension of the three dimensional array
3256: . mstart - first index you will use in first coordinate direction (often 0)
3257: . nstart - first index in the second coordinate direction (often 0)
3258: . pstart - first index in the third coordinate direction (often 0)
3259: - a - location of pointer to array obtained from VecGetArray3d()
3261: Level: developer
3263: Notes:
3264: For regular PETSc vectors this routine does not involve any copies. For
3265: any special vectors that do not store local vector data in a contiguous
3266: array, this routine will copy the data back into the underlying
3267: vector data structure from the array obtained with `VecGetArray()`.
3269: This routine actually zeros out the `a` pointer.
3271: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3272: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3273: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3274: @*/
3275: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3276: {
3277: void *dummy;
3279: PetscFunctionBegin;
3281: PetscAssertPointer(a, 8);
3283: dummy = (void *)(*a + mstart);
3284: PetscCall(PetscFree(dummy));
3285: PetscCall(VecRestoreArrayWrite(x, NULL));
3286: *a = NULL;
3287: PetscFunctionReturn(PETSC_SUCCESS);
3288: }
3290: /*@C
3291: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3292: processor's portion of the vector data. You MUST call `VecRestoreArray4d()`
3293: when you no longer need access to the array.
3295: Logically Collective
3297: Input Parameters:
3298: + x - the vector
3299: . m - first dimension of four dimensional array
3300: . n - second dimension of four dimensional array
3301: . p - third dimension of four dimensional array
3302: . q - fourth dimension of four dimensional array
3303: . mstart - first index you will use in first coordinate direction (often 0)
3304: . nstart - first index in the second coordinate direction (often 0)
3305: . pstart - first index in the third coordinate direction (often 0)
3306: - qstart - first index in the fourth coordinate direction (often 0)
3308: Output Parameter:
3309: . a - location to put pointer to the array
3311: Level: developer
3313: Notes:
3314: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3315: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3316: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3317: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3319: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3321: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3322: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3323: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3324: @*/
3325: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3326: {
3327: PetscInt i, N, j, k;
3328: PetscScalar *aa, ***b, **c;
3330: PetscFunctionBegin;
3332: PetscAssertPointer(a, 10);
3334: PetscCall(VecGetLocalSize(x, &N));
3335: 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);
3336: PetscCall(VecGetArray(x, &aa));
3338: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3339: b = (PetscScalar ***)((*a) + m);
3340: c = (PetscScalar **)(b + m * n);
3341: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3342: for (i = 0; i < m; i++)
3343: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3344: for (i = 0; i < m; i++)
3345: for (j = 0; j < n; j++)
3346: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3347: *a -= mstart;
3348: PetscFunctionReturn(PETSC_SUCCESS);
3349: }
3351: /*@C
3352: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3353: processor's portion of the vector data. You MUST call `VecRestoreArray4dWrite()`
3354: when you no longer need access to the array.
3356: Logically Collective
3358: Input Parameters:
3359: + x - the vector
3360: . m - first dimension of four dimensional array
3361: . n - second dimension of four dimensional array
3362: . p - third dimension of four dimensional array
3363: . q - fourth dimension of four dimensional array
3364: . mstart - first index you will use in first coordinate direction (often 0)
3365: . nstart - first index in the second coordinate direction (often 0)
3366: . pstart - first index in the third coordinate direction (often 0)
3367: - qstart - first index in the fourth coordinate direction (often 0)
3369: Output Parameter:
3370: . a - location to put pointer to the array
3372: Level: developer
3374: Notes:
3375: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3376: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3377: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3378: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3380: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3382: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3383: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3384: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3385: @*/
3386: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3387: {
3388: PetscInt i, N, j, k;
3389: PetscScalar *aa, ***b, **c;
3391: PetscFunctionBegin;
3393: PetscAssertPointer(a, 10);
3395: PetscCall(VecGetLocalSize(x, &N));
3396: 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);
3397: PetscCall(VecGetArrayWrite(x, &aa));
3399: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3400: b = (PetscScalar ***)((*a) + m);
3401: c = (PetscScalar **)(b + m * n);
3402: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3403: for (i = 0; i < m; i++)
3404: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3405: for (i = 0; i < m; i++)
3406: for (j = 0; j < n; j++)
3407: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3408: *a -= mstart;
3409: PetscFunctionReturn(PETSC_SUCCESS);
3410: }
3412: /*@C
3413: VecRestoreArray4d - Restores a vector after `VecGetArray4d()` has been called.
3415: Logically Collective
3417: Input Parameters:
3418: + x - the vector
3419: . m - first dimension of four dimensional array
3420: . n - second dimension of the four dimensional array
3421: . p - third dimension of the four dimensional array
3422: . q - fourth dimension of the four dimensional array
3423: . mstart - first index you will use in first coordinate direction (often 0)
3424: . nstart - first index in the second coordinate direction (often 0)
3425: . pstart - first index in the third coordinate direction (often 0)
3426: . qstart - first index in the fourth coordinate direction (often 0)
3427: - a - location of pointer to array obtained from VecGetArray4d()
3429: Level: developer
3431: Notes:
3432: For regular PETSc vectors this routine does not involve any copies. For
3433: any special vectors that do not store local vector data in a contiguous
3434: array, this routine will copy the data back into the underlying
3435: vector data structure from the array obtained with `VecGetArray()`.
3437: This routine actually zeros out the `a` pointer.
3439: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3440: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3441: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecGet`
3442: @*/
3443: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3444: {
3445: void *dummy;
3447: PetscFunctionBegin;
3449: PetscAssertPointer(a, 10);
3451: dummy = (void *)(*a + mstart);
3452: PetscCall(PetscFree(dummy));
3453: PetscCall(VecRestoreArray(x, NULL));
3454: *a = NULL;
3455: PetscFunctionReturn(PETSC_SUCCESS);
3456: }
3458: /*@C
3459: VecRestoreArray4dWrite - Restores a vector after `VecGetArray4dWrite()` has been called.
3461: Logically Collective
3463: Input Parameters:
3464: + x - the vector
3465: . m - first dimension of four dimensional array
3466: . n - second dimension of the four dimensional array
3467: . p - third dimension of the four dimensional array
3468: . q - fourth dimension of the four dimensional array
3469: . mstart - first index you will use in first coordinate direction (often 0)
3470: . nstart - first index in the second coordinate direction (often 0)
3471: . pstart - first index in the third coordinate direction (often 0)
3472: . qstart - first index in the fourth coordinate direction (often 0)
3473: - a - location of pointer to array obtained from `VecGetArray4d()`
3475: Level: developer
3477: Notes:
3478: For regular PETSc vectors this routine does not involve any copies. For
3479: any special vectors that do not store local vector data in a contiguous
3480: array, this routine will copy the data back into the underlying
3481: vector data structure from the array obtained with `VecGetArray()`.
3483: This routine actually zeros out the `a` pointer.
3485: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3486: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3487: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3488: @*/
3489: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3490: {
3491: void *dummy;
3493: PetscFunctionBegin;
3495: PetscAssertPointer(a, 10);
3497: dummy = (void *)(*a + mstart);
3498: PetscCall(PetscFree(dummy));
3499: PetscCall(VecRestoreArrayWrite(x, NULL));
3500: *a = NULL;
3501: PetscFunctionReturn(PETSC_SUCCESS);
3502: }
3504: /*@C
3505: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3506: processor's portion of the vector data. You MUST call `VecRestoreArray2dRead()`
3507: when you no longer need access to the array.
3509: Logically Collective
3511: Input Parameters:
3512: + x - the vector
3513: . m - first dimension of two dimensional array
3514: . n - second dimension of two dimensional array
3515: . mstart - first index you will use in first coordinate direction (often 0)
3516: - nstart - first index in the second coordinate direction (often 0)
3518: Output Parameter:
3519: . a - location to put pointer to the array
3521: Level: developer
3523: Notes:
3524: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
3525: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3526: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3527: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
3529: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3531: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3532: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3533: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3534: @*/
3535: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3536: {
3537: PetscInt i, N;
3538: const PetscScalar *aa;
3540: PetscFunctionBegin;
3542: PetscAssertPointer(a, 6);
3544: PetscCall(VecGetLocalSize(x, &N));
3545: 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);
3546: PetscCall(VecGetArrayRead(x, &aa));
3548: PetscCall(PetscMalloc1(m, a));
3549: for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3550: *a -= mstart;
3551: PetscFunctionReturn(PETSC_SUCCESS);
3552: }
3554: /*@C
3555: VecRestoreArray2dRead - Restores a vector after `VecGetArray2dRead()` has been called.
3557: Logically Collective
3559: Input Parameters:
3560: + x - the vector
3561: . m - first dimension of two dimensional array
3562: . n - second dimension of the two dimensional array
3563: . mstart - first index you will use in first coordinate direction (often 0)
3564: . nstart - first index in the second coordinate direction (often 0)
3565: - a - location of pointer to array obtained from VecGetArray2d()
3567: Level: developer
3569: Notes:
3570: For regular PETSc vectors this routine does not involve any copies. For
3571: any special vectors that do not store local vector data in a contiguous
3572: array, this routine will copy the data back into the underlying
3573: vector data structure from the array obtained with `VecGetArray()`.
3575: This routine actually zeros out the `a` pointer.
3577: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3578: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3579: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3580: @*/
3581: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3582: {
3583: void *dummy;
3585: PetscFunctionBegin;
3587: PetscAssertPointer(a, 6);
3589: dummy = (void *)(*a + mstart);
3590: PetscCall(PetscFree(dummy));
3591: PetscCall(VecRestoreArrayRead(x, NULL));
3592: *a = NULL;
3593: PetscFunctionReturn(PETSC_SUCCESS);
3594: }
3596: /*@C
3597: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3598: processor's portion of the vector data. You MUST call `VecRestoreArray1dRead()`
3599: when you no longer need access to the array.
3601: Logically Collective
3603: Input Parameters:
3604: + x - the vector
3605: . m - first dimension of two dimensional array
3606: - mstart - first index you will use in first coordinate direction (often 0)
3608: Output Parameter:
3609: . a - location to put pointer to the array
3611: Level: developer
3613: Notes:
3614: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3615: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3616: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3618: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3620: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3621: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3622: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3623: @*/
3624: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3625: {
3626: PetscInt N;
3628: PetscFunctionBegin;
3630: PetscAssertPointer(a, 4);
3632: PetscCall(VecGetLocalSize(x, &N));
3633: 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);
3634: PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3635: *a -= mstart;
3636: PetscFunctionReturn(PETSC_SUCCESS);
3637: }
3639: /*@C
3640: VecRestoreArray1dRead - Restores a vector after `VecGetArray1dRead()` has been called.
3642: Logically Collective
3644: Input Parameters:
3645: + x - the vector
3646: . m - first dimension of two dimensional array
3647: . mstart - first index you will use in first coordinate direction (often 0)
3648: - a - location of pointer to array obtained from `VecGetArray1dRead()`
3650: Level: developer
3652: Notes:
3653: For regular PETSc vectors this routine does not involve any copies. For
3654: any special vectors that do not store local vector data in a contiguous
3655: array, this routine will copy the data back into the underlying
3656: vector data structure from the array obtained with `VecGetArray1dRead()`.
3658: This routine actually zeros out the `a` pointer.
3660: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3661: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3662: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3663: @*/
3664: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3665: {
3666: PetscFunctionBegin;
3669: PetscCall(VecRestoreArrayRead(x, NULL));
3670: *a = NULL;
3671: PetscFunctionReturn(PETSC_SUCCESS);
3672: }
3674: /*@C
3675: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3676: processor's portion of the vector data. You MUST call `VecRestoreArray3dRead()`
3677: when you no longer need access to the array.
3679: Logically Collective
3681: Input Parameters:
3682: + x - the vector
3683: . m - first dimension of three dimensional array
3684: . n - second dimension of three dimensional array
3685: . p - third dimension of three dimensional array
3686: . mstart - first index you will use in first coordinate direction (often 0)
3687: . nstart - first index in the second coordinate direction (often 0)
3688: - pstart - first index in the third coordinate direction (often 0)
3690: Output Parameter:
3691: . a - location to put pointer to the array
3693: Level: developer
3695: Notes:
3696: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3697: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3698: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3699: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3dRead()`.
3701: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3703: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3704: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3705: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3706: @*/
3707: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3708: {
3709: PetscInt i, N, j;
3710: const PetscScalar *aa;
3711: PetscScalar **b;
3713: PetscFunctionBegin;
3715: PetscAssertPointer(a, 8);
3717: PetscCall(VecGetLocalSize(x, &N));
3718: 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);
3719: PetscCall(VecGetArrayRead(x, &aa));
3721: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3722: b = (PetscScalar **)((*a) + m);
3723: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3724: for (i = 0; i < m; i++)
3725: for (j = 0; j < n; j++) b[i * n + j] = (PetscScalar *)aa + i * n * p + j * p - pstart;
3726: *a -= mstart;
3727: PetscFunctionReturn(PETSC_SUCCESS);
3728: }
3730: /*@C
3731: VecRestoreArray3dRead - Restores a vector after `VecGetArray3dRead()` has been called.
3733: Logically Collective
3735: Input Parameters:
3736: + x - the vector
3737: . m - first dimension of three dimensional array
3738: . n - second dimension of the three dimensional array
3739: . p - third dimension of the three dimensional array
3740: . mstart - first index you will use in first coordinate direction (often 0)
3741: . nstart - first index in the second coordinate direction (often 0)
3742: . pstart - first index in the third coordinate direction (often 0)
3743: - a - location of pointer to array obtained from `VecGetArray3dRead()`
3745: Level: developer
3747: Notes:
3748: For regular PETSc vectors this routine does not involve any copies. For
3749: any special vectors that do not store local vector data in a contiguous
3750: array, this routine will copy the data back into the underlying
3751: vector data structure from the array obtained with `VecGetArray()`.
3753: This routine actually zeros out the `a` pointer.
3755: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3756: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3757: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3758: @*/
3759: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3760: {
3761: void *dummy;
3763: PetscFunctionBegin;
3765: PetscAssertPointer(a, 8);
3767: dummy = (void *)(*a + mstart);
3768: PetscCall(PetscFree(dummy));
3769: PetscCall(VecRestoreArrayRead(x, NULL));
3770: *a = NULL;
3771: PetscFunctionReturn(PETSC_SUCCESS);
3772: }
3774: /*@C
3775: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3776: processor's portion of the vector data. You MUST call `VecRestoreArray4dRead()`
3777: when you no longer need access to the array.
3779: Logically Collective
3781: Input Parameters:
3782: + x - the vector
3783: . m - first dimension of four dimensional array
3784: . n - second dimension of four dimensional array
3785: . p - third dimension of four dimensional array
3786: . q - fourth dimension of four dimensional array
3787: . mstart - first index you will use in first coordinate direction (often 0)
3788: . nstart - first index in the second coordinate direction (often 0)
3789: . pstart - first index in the third coordinate direction (often 0)
3790: - qstart - first index in the fourth coordinate direction (often 0)
3792: Output Parameter:
3793: . a - location to put pointer to the array
3795: Level: beginner
3797: Notes:
3798: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3799: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3800: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3801: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3803: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3805: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3806: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3807: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3808: @*/
3809: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3810: {
3811: PetscInt i, N, j, k;
3812: const PetscScalar *aa;
3813: PetscScalar ***b, **c;
3815: PetscFunctionBegin;
3817: PetscAssertPointer(a, 10);
3819: PetscCall(VecGetLocalSize(x, &N));
3820: 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);
3821: PetscCall(VecGetArrayRead(x, &aa));
3823: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3824: b = (PetscScalar ***)((*a) + m);
3825: c = (PetscScalar **)(b + m * n);
3826: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3827: for (i = 0; i < m; i++)
3828: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3829: for (i = 0; i < m; i++)
3830: for (j = 0; j < n; j++)
3831: 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;
3832: *a -= mstart;
3833: PetscFunctionReturn(PETSC_SUCCESS);
3834: }
3836: /*@C
3837: VecRestoreArray4dRead - Restores a vector after `VecGetArray4d()` has been called.
3839: Logically Collective
3841: Input Parameters:
3842: + x - the vector
3843: . m - first dimension of four dimensional array
3844: . n - second dimension of the four dimensional array
3845: . p - third dimension of the four dimensional array
3846: . q - fourth dimension of the four dimensional array
3847: . mstart - first index you will use in first coordinate direction (often 0)
3848: . nstart - first index in the second coordinate direction (often 0)
3849: . pstart - first index in the third coordinate direction (often 0)
3850: . qstart - first index in the fourth coordinate direction (often 0)
3851: - a - location of pointer to array obtained from `VecGetArray4dRead()`
3853: Level: beginner
3855: Notes:
3856: For regular PETSc vectors this routine does not involve any copies. For
3857: any special vectors that do not store local vector data in a contiguous
3858: array, this routine will copy the data back into the underlying
3859: vector data structure from the array obtained with `VecGetArray()`.
3861: This routine actually zeros out the `a` pointer.
3863: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3864: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3865: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3866: @*/
3867: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3868: {
3869: void *dummy;
3871: PetscFunctionBegin;
3873: PetscAssertPointer(a, 10);
3875: dummy = (void *)(*a + mstart);
3876: PetscCall(PetscFree(dummy));
3877: PetscCall(VecRestoreArrayRead(x, NULL));
3878: *a = NULL;
3879: PetscFunctionReturn(PETSC_SUCCESS);
3880: }
3882: #if defined(PETSC_USE_DEBUG)
3884: /*@
3885: VecLockGet - Gets the current lock status of a vector
3887: Logically Collective
3889: Input Parameter:
3890: . x - the vector
3892: Output Parameter:
3893: . state - greater than zero indicates the vector is locked for read; less than zero indicates the vector is
3894: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3896: Level: advanced
3898: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3899: @*/
3900: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3901: {
3902: PetscFunctionBegin;
3904: *state = x->lock;
3905: PetscFunctionReturn(PETSC_SUCCESS);
3906: }
3908: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3909: {
3910: PetscFunctionBegin;
3912: PetscAssertPointer(file, 2);
3913: PetscAssertPointer(func, 3);
3914: PetscAssertPointer(line, 4);
3915: #if !PetscDefined(HAVE_THREADSAFETY)
3916: {
3917: const int index = x->lockstack.currentsize - 1;
3919: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupted vec lock stack, have negative index %d", index);
3920: *file = x->lockstack.file[index];
3921: *func = x->lockstack.function[index];
3922: *line = x->lockstack.line[index];
3923: }
3924: #else
3925: *file = NULL;
3926: *func = NULL;
3927: *line = 0;
3928: #endif
3929: PetscFunctionReturn(PETSC_SUCCESS);
3930: }
3932: /*@
3933: VecLockReadPush - Pushes a read-only lock on a vector to prevent it from being written to
3935: Logically Collective
3937: Input Parameter:
3938: . x - the vector
3940: Level: intermediate
3942: Notes:
3943: If this is set then calls to `VecGetArray()` or `VecSetValues()` or any other routines that change the vectors values will generate an error.
3945: The call can be nested, i.e., called multiple times on the same vector, but each `VecLockReadPush()` has to have one matching
3946: `VecLockReadPop()`, which removes the latest read-only lock.
3948: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3949: @*/
3950: PetscErrorCode VecLockReadPush(Vec x)
3951: {
3952: PetscFunctionBegin;
3954: 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");
3955: #if !PetscDefined(HAVE_THREADSAFETY)
3956: {
3957: const char *file, *func;
3958: int index, line;
3960: if ((index = petscstack.currentsize - 2) == -1) {
3961: // vec was locked "outside" of petsc, either in user-land or main. the error message will
3962: // now show this function as the culprit, but it will include the stacktrace
3963: file = "unknown user-file";
3964: func = "unknown_user_function";
3965: line = 0;
3966: } else {
3967: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected petscstack, have negative index %d", index);
3968: file = petscstack.file[index];
3969: func = petscstack.function[index];
3970: line = petscstack.line[index];
3971: }
3972: PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
3973: }
3974: #endif
3975: PetscFunctionReturn(PETSC_SUCCESS);
3976: }
3978: /*@
3979: VecLockReadPop - Pops a read-only lock from a vector
3981: Logically Collective
3983: Input Parameter:
3984: . x - the vector
3986: Level: intermediate
3988: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
3989: @*/
3990: PetscErrorCode VecLockReadPop(Vec x)
3991: {
3992: PetscFunctionBegin;
3994: PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
3995: #if !PetscDefined(HAVE_THREADSAFETY)
3996: {
3997: const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];
3999: PetscStackPop_Private(x->lockstack, previous);
4000: }
4001: #endif
4002: PetscFunctionReturn(PETSC_SUCCESS);
4003: }
4005: /*@C
4006: VecLockWriteSet - Lock or unlock a vector for exclusive read/write access
4008: Logically Collective
4010: Input Parameters:
4011: + x - the vector
4012: - flg - `PETSC_TRUE` to lock the vector for exclusive read/write access; `PETSC_FALSE` to unlock it.
4014: Level: intermediate
4016: Notes:
4017: The function is useful in split-phase computations, which usually have a begin phase and an end phase.
4018: One can call `VecLockWriteSet`(x,`PETSC_TRUE`) in the begin phase to lock a vector for exclusive
4019: access, and call `VecLockWriteSet`(x,`PETSC_FALSE`) in the end phase to unlock the vector from exclusive
4020: access. In this way, one is ensured no other operations can access the vector in between. The code may like
4022: .vb
4023: VecGetArray(x,&xdata); // begin phase
4024: VecLockWriteSet(v,PETSC_TRUE);
4026: Other operations, which can not access x anymore (they can access xdata, of course)
4028: VecRestoreArray(x,&vdata); // end phase
4029: VecLockWriteSet(v,PETSC_FALSE);
4030: .ve
4032: The call can not be nested on the same vector, in other words, one can not call `VecLockWriteSet`(x,`PETSC_TRUE`)
4033: again before calling `VecLockWriteSet`(v,`PETSC_FALSE`).
4035: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
4036: @*/
4037: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
4038: {
4039: PetscFunctionBegin;
4041: if (flg) {
4042: 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");
4043: 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");
4044: x->lock = -1;
4045: } else {
4046: 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");
4047: x->lock = 0;
4048: }
4049: PetscFunctionReturn(PETSC_SUCCESS);
4050: }
4052: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
4053: /*@
4054: VecLockPush - Pushes a read-only lock on a vector to prevent it from being written to
4056: Level: deprecated
4058: .seealso: [](ch_vectors), `Vec`, `VecLockReadPush()`
4059: @*/
4060: PetscErrorCode VecLockPush(Vec x)
4061: {
4062: PetscFunctionBegin;
4063: PetscCall(VecLockReadPush(x));
4064: PetscFunctionReturn(PETSC_SUCCESS);
4065: }
4067: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
4068: /*@
4069: VecLockPop - Pops a read-only lock from a vector
4071: Level: deprecated
4073: .seealso: [](ch_vectors), `Vec`, `VecLockReadPop()`
4074: @*/
4075: PetscErrorCode VecLockPop(Vec x)
4076: {
4077: PetscFunctionBegin;
4078: PetscCall(VecLockReadPop(x));
4079: PetscFunctionReturn(PETSC_SUCCESS);
4080: }
4082: #endif