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: /* Cached data? */
209: PetscCall(VecNormAvailable(x, type, &flg, val));
210: if (flg) PetscFunctionReturn(PETSC_SUCCESS);
212: PetscCall(VecLockReadPush(x));
213: PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
214: PetscUseTypeMethod(x, norm, type, val);
215: PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
216: PetscCall(VecLockReadPop(x));
218: if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: /*@
223: VecNormAvailable - Returns the vector norm if it is already known. That is, it has been previously computed and cached in the vector
225: Not Collective
227: Input Parameters:
228: + x - the vector
229: - 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
230: `NORM_1_AND_2`, which computes both norms and stores them
231: in a two element array.
233: Output Parameters:
234: + available - `PETSC_TRUE` if the val returned is valid
235: - val - the norm
237: Level: intermediate
239: Developer Notes:
240: `PETSC_HAVE_SLOW_BLAS_NORM2` will cause a C (loop unrolled) version of the norm to be used, rather
241: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
242: (as opposed to vendor provided) because the FORTRAN BLAS `NRM2()` routine is very slow.
244: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
245: `VecNormBegin()`, `VecNormEnd()`
246: @*/
247: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
248: {
249: PetscFunctionBegin;
252: PetscAssertPointer(available, 3);
253: PetscAssertPointer(val, 4);
255: if (type == NORM_1_AND_2) {
256: *available = PETSC_FALSE;
257: } else {
258: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
259: }
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: /*@
264: VecNormalize - Normalizes a vector by its 2-norm.
266: Collective
268: Input Parameter:
269: . x - the vector
271: Output Parameter:
272: . val - the vector norm before normalization. May be `NULL` if the value is not needed.
274: Level: intermediate
276: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
277: @*/
278: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
279: {
280: PetscReal norm;
282: PetscFunctionBegin;
285: PetscCall(VecSetErrorIfLocked(x, 1));
286: if (val) PetscAssertPointer(val, 2);
287: PetscCall(PetscLogEventBegin(VEC_Normalize, x, 0, 0, 0));
288: PetscCall(VecNorm(x, NORM_2, &norm));
289: if (norm == 0.0) PetscCall(PetscInfo(x, "Vector of zero norm can not be normalized; Returning only the zero norm\n"));
290: else if (PetscIsInfOrNanReal(norm)) PetscCall(PetscInfo(x, "Vector with Inf or Nan norm can not be normalized; Returning only the norm\n"));
291: else {
292: PetscScalar s = 1.0 / norm;
293: PetscCall(VecScale(x, s));
294: }
295: PetscCall(PetscLogEventEnd(VEC_Normalize, x, 0, 0, 0));
296: if (val) *val = norm;
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: /*@C
301: VecMax - Determines the vector component with maximum real part and its location.
303: Collective
305: Input Parameter:
306: . x - the vector
308: Output Parameters:
309: + p - the index of `val` (pass `NULL` if you don't want this) in the vector
310: - val - the maximum component
312: Level: intermediate
314: Notes:
315: Returns the value `PETSC_MIN_REAL` and negative `p` if the vector is of length 0.
317: Returns the smallest index with the maximum value
319: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
320: @*/
321: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
322: {
323: PetscFunctionBegin;
326: VecCheckAssembled(x);
327: if (p) PetscAssertPointer(p, 2);
328: PetscAssertPointer(val, 3);
329: PetscCall(VecLockReadPush(x));
330: PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
331: PetscUseTypeMethod(x, max, p, val);
332: PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
333: PetscCall(VecLockReadPop(x));
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: /*@C
338: VecMin - Determines the vector component with minimum real part and its location.
340: Collective
342: Input Parameter:
343: . x - the vector
345: Output Parameters:
346: + p - the index of `val` (pass `NULL` if you don't want this location) in the vector
347: - val - the minimum component
349: Level: intermediate
351: Notes:
352: Returns the value `PETSC_MAX_REAL` and negative `p` if the vector is of length 0.
354: This returns the smallest index with the minimum value
356: .seealso: [](ch_vectors), `Vec`, `VecMax()`
357: @*/
358: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
359: {
360: PetscFunctionBegin;
363: VecCheckAssembled(x);
364: if (p) PetscAssertPointer(p, 2);
365: PetscAssertPointer(val, 3);
366: PetscCall(VecLockReadPush(x));
367: PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
368: PetscUseTypeMethod(x, min, p, val);
369: PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
370: PetscCall(VecLockReadPop(x));
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: /*@
375: VecTDot - Computes an indefinite vector dot product. That is, this
376: routine does NOT use the complex conjugate.
378: Collective
380: Input Parameters:
381: + x - first vector
382: - y - second vector
384: Output Parameter:
385: . val - the dot product
387: Level: intermediate
389: Notes for Users of Complex Numbers:
390: For complex vectors, `VecTDot()` computes the indefinite form
391: $ val = (x,y) = y^T x,
392: where y^T denotes the transpose of y.
394: Use `VecDot()` for the inner product
395: $ val = (x,y) = y^H x,
396: where y^H denotes the conjugate transpose of y.
398: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
399: @*/
400: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
401: {
402: PetscFunctionBegin;
405: PetscAssertPointer(val, 3);
408: PetscCheckSameTypeAndComm(x, 1, y, 2);
409: VecCheckSameSize(x, 1, y, 2);
410: VecCheckAssembled(x);
411: VecCheckAssembled(y);
413: PetscCall(VecLockReadPush(x));
414: PetscCall(VecLockReadPush(y));
415: PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
416: PetscUseTypeMethod(x, tdot, y, val);
417: PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
418: PetscCall(VecLockReadPop(x));
419: PetscCall(VecLockReadPop(y));
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: PetscErrorCode VecScaleAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
424: {
425: PetscReal norms[4];
426: PetscBool flgs[4];
427: PetscScalar one = 1.0;
429: PetscFunctionBegin;
432: VecCheckAssembled(x);
433: PetscCall(VecSetErrorIfLocked(x, 1));
434: if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);
436: /* get current stashed norms */
437: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
439: PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
440: VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
441: PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));
443: PetscCall(PetscObjectStateIncrease((PetscObject)x));
444: /* put the scaled stashed norms back into the Vec */
445: for (PetscInt i = 0; i < 4; i++) {
446: PetscReal ar = PetscAbsScalar(alpha);
447: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
448: }
449: PetscFunctionReturn(PETSC_SUCCESS);
450: }
452: /*@
453: VecScale - Scales a vector.
455: Not Collective
457: Input Parameters:
458: + x - the vector
459: - alpha - the scalar
461: Level: intermediate
463: Note:
464: For a vector with n components, `VecScale()` computes x[i] = alpha * x[i], for i=1,...,n.
466: .seealso: [](ch_vectors), `Vec`, `VecSet()`
467: @*/
468: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
469: {
470: PetscFunctionBegin;
471: PetscCall(VecScaleAsync_Private(x, alpha, NULL));
472: PetscFunctionReturn(PETSC_SUCCESS);
473: }
475: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
476: {
477: PetscFunctionBegin;
480: VecCheckAssembled(x);
482: PetscCall(VecSetErrorIfLocked(x, 1));
484: if (alpha == 0) {
485: PetscReal norm;
486: PetscBool set;
488: PetscCall(VecNormAvailable(x, NORM_2, &set, &norm));
489: if (set == PETSC_TRUE && norm == 0) PetscFunctionReturn(PETSC_SUCCESS);
490: }
491: PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
492: VecMethodDispatch(x, dctx, VecAsyncFnName(Set), set, (Vec, PetscScalar, PetscDeviceContext), alpha);
493: PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
494: PetscCall(PetscObjectStateIncrease((PetscObject)x));
496: /* norms can be simply set (if |alpha|*N not too large) */
498: {
499: PetscReal val = PetscAbsScalar(alpha);
500: const PetscInt N = x->map->N;
502: if (N == 0) {
503: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0l));
504: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
505: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
506: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
507: } else if (val > PETSC_MAX_REAL / N) {
508: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
509: } else {
510: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
511: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
512: val *= PetscSqrtReal((PetscReal)N);
513: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
514: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
515: }
516: }
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }
520: /*@
521: VecSet - Sets all components of a vector to a single scalar value.
523: Logically Collective
525: Input Parameters:
526: + x - the vector
527: - alpha - the scalar
529: Level: beginner
531: Notes:
532: For a vector of dimension n, `VecSet()` sets x[i] = alpha, for i=1,...,n,
533: so that all vector entries then equal the identical
534: scalar value, `alpha`. Use the more general routine
535: `VecSetValues()` to set different vector entries.
537: You CANNOT call this after you have called `VecSetValues()` but before you call
538: `VecAssemblyBegin()`
540: If `alpha` is zero and the norm of the vector is known to be zero then this skips the unneeded zeroing process
542: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
543: @*/
544: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
545: {
546: PetscFunctionBegin;
547: PetscCall(VecSetAsync_Private(x, alpha, NULL));
548: PetscFunctionReturn(PETSC_SUCCESS);
549: }
551: PetscErrorCode VecAXPYAsync_Private(Vec y, PetscScalar alpha, Vec x, PetscDeviceContext dctx)
552: {
553: PetscFunctionBegin;
558: PetscCheckSameTypeAndComm(x, 3, y, 1);
559: VecCheckSameSize(x, 3, y, 1);
560: VecCheckAssembled(x);
561: VecCheckAssembled(y);
563: if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
564: PetscCall(VecSetErrorIfLocked(y, 1));
565: if (x == y) {
566: PetscCall(VecScale(y, alpha + 1.0));
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
569: PetscCall(VecLockReadPush(x));
570: PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
571: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPY), axpy, (Vec, PetscScalar, Vec, PetscDeviceContext), alpha, x);
572: PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
573: PetscCall(VecLockReadPop(x));
574: PetscCall(PetscObjectStateIncrease((PetscObject)y));
575: PetscFunctionReturn(PETSC_SUCCESS);
576: }
577: /*@
578: VecAXPY - Computes `y = alpha x + y`.
580: Logically Collective
582: Input Parameters:
583: + alpha - the scalar
584: . x - vector scale by `alpha`
585: - y - vector accumulated into
587: Output Parameter:
588: . y - output vector
590: Level: intermediate
592: Notes:
593: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
594: .vb
595: VecAXPY(y,alpha,x) y = alpha x + y
596: VecAYPX(y,beta,x) y = x + beta y
597: VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
598: VecWAXPY(w,alpha,x,y) w = alpha x + y
599: VecAXPBYPCZ(z,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
600: VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
601: .ve
603: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
604: @*/
605: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
606: {
607: PetscFunctionBegin;
608: PetscCall(VecAXPYAsync_Private(y, alpha, x, NULL));
609: PetscFunctionReturn(PETSC_SUCCESS);
610: }
612: PetscErrorCode VecAYPXAsync_Private(Vec y, PetscScalar beta, Vec x, PetscDeviceContext dctx)
613: {
614: PetscFunctionBegin;
619: PetscCheckSameTypeAndComm(x, 3, y, 1);
620: VecCheckSameSize(x, 1, y, 3);
621: VecCheckAssembled(x);
622: VecCheckAssembled(y);
624: PetscCall(VecSetErrorIfLocked(y, 1));
625: if (x == y) {
626: PetscCall(VecScale(y, beta + 1.0));
627: PetscFunctionReturn(PETSC_SUCCESS);
628: }
629: PetscCall(VecLockReadPush(x));
630: if (beta == (PetscScalar)0.0) {
631: PetscCall(VecCopy(x, y));
632: } else {
633: PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
634: VecMethodDispatch(y, dctx, VecAsyncFnName(AYPX), aypx, (Vec, PetscScalar, Vec, PetscDeviceContext), beta, x);
635: PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
636: PetscCall(PetscObjectStateIncrease((PetscObject)y));
637: }
638: PetscCall(VecLockReadPop(x));
639: PetscFunctionReturn(PETSC_SUCCESS);
640: }
642: /*@
643: VecAYPX - Computes `y = x + beta y`.
645: Logically Collective
647: Input Parameters:
648: + beta - the scalar
649: . x - the unscaled vector
650: - y - the vector to be scaled
652: Output Parameter:
653: . y - output vector
655: Level: intermediate
657: Developer Notes:
658: The implementation is optimized for `beta` of -1.0, 0.0, and 1.0
660: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
661: @*/
662: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
663: {
664: PetscFunctionBegin;
665: PetscCall(VecAYPXAsync_Private(y, beta, x, NULL));
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: PetscErrorCode VecAXPBYAsync_Private(Vec y, PetscScalar alpha, PetscScalar beta, Vec x, PetscDeviceContext dctx)
670: {
671: PetscFunctionBegin;
676: PetscCheckSameTypeAndComm(x, 4, y, 1);
677: VecCheckSameSize(y, 1, x, 4);
678: VecCheckAssembled(x);
679: VecCheckAssembled(y);
682: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
683: if (x == y) {
684: PetscCall(VecScale(y, alpha + beta));
685: PetscFunctionReturn(PETSC_SUCCESS);
686: }
688: PetscCall(VecSetErrorIfLocked(y, 1));
689: PetscCall(VecLockReadPush(x));
690: PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
691: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPBY), axpby, (Vec, PetscScalar, PetscScalar, Vec, PetscDeviceContext), alpha, beta, x);
692: PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
693: PetscCall(PetscObjectStateIncrease((PetscObject)y));
694: PetscCall(VecLockReadPop(x));
695: PetscFunctionReturn(PETSC_SUCCESS);
696: }
697: /*@
698: VecAXPBY - Computes `y = alpha x + beta y`.
700: Logically Collective
702: Input Parameters:
703: + alpha - first scalar
704: . beta - second scalar
705: . x - the first scaled vector
706: - y - the second scaled vector
708: Output Parameter:
709: . y - output vector
711: Level: intermediate
713: Developer Notes:
714: The implementation is optimized for `alpha` and/or `beta` values of 0.0 and 1.0
716: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
717: @*/
718: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
719: {
720: PetscFunctionBegin;
721: PetscCall(VecAXPBYAsync_Private(y, alpha, beta, x, NULL));
722: PetscFunctionReturn(PETSC_SUCCESS);
723: }
725: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
726: {
727: PetscFunctionBegin;
734: PetscCheckSameTypeAndComm(x, 5, y, 6);
735: PetscCheckSameTypeAndComm(x, 5, z, 1);
736: VecCheckSameSize(x, 1, y, 5);
737: VecCheckSameSize(x, 1, z, 6);
738: PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
739: PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
740: VecCheckAssembled(x);
741: VecCheckAssembled(y);
742: VecCheckAssembled(z);
746: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
748: PetscCall(VecSetErrorIfLocked(z, 1));
749: PetscCall(VecLockReadPush(x));
750: PetscCall(VecLockReadPush(y));
751: PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
752: VecMethodDispatch(z, dctx, VecAsyncFnName(AXPBYPCZ), axpbypcz, (Vec, PetscScalar, PetscScalar, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, beta, gamma, x, y);
753: PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
754: PetscCall(PetscObjectStateIncrease((PetscObject)z));
755: PetscCall(VecLockReadPop(x));
756: PetscCall(VecLockReadPop(y));
757: PetscFunctionReturn(PETSC_SUCCESS);
758: }
759: /*@
760: VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`
762: Logically Collective
764: Input Parameters:
765: + alpha - first scalar
766: . beta - second scalar
767: . gamma - third scalar
768: . x - first vector
769: . y - second vector
770: - z - third vector
772: Output Parameter:
773: . z - output vector
775: Level: intermediate
777: Note:
778: `x`, `y` and `z` must be different vectors
780: Developer Notes:
781: The implementation is optimized for `alpha` of 1.0 and `gamma` of 1.0 or 0.0
783: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
784: @*/
785: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
786: {
787: PetscFunctionBegin;
788: PetscCall(VecAXPBYPCZAsync_Private(z, alpha, beta, gamma, x, y, NULL));
789: PetscFunctionReturn(PETSC_SUCCESS);
790: }
792: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
793: {
794: PetscFunctionBegin;
801: PetscCheckSameTypeAndComm(x, 3, y, 4);
802: PetscCheckSameTypeAndComm(y, 4, w, 1);
803: VecCheckSameSize(x, 3, y, 4);
804: VecCheckSameSize(x, 3, w, 1);
805: PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
806: PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
807: VecCheckAssembled(x);
808: VecCheckAssembled(y);
810: PetscCall(VecSetErrorIfLocked(w, 1));
812: PetscCall(VecLockReadPush(x));
813: PetscCall(VecLockReadPush(y));
814: if (alpha == (PetscScalar)0.0) {
815: PetscCall(VecCopyAsync_Private(y, w, dctx));
816: } else {
817: PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
818: VecMethodDispatch(w, dctx, VecAsyncFnName(WAXPY), waxpy, (Vec, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, x, y);
819: PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
820: PetscCall(PetscObjectStateIncrease((PetscObject)w));
821: }
822: PetscCall(VecLockReadPop(x));
823: PetscCall(VecLockReadPop(y));
824: PetscFunctionReturn(PETSC_SUCCESS);
825: }
827: /*@
828: VecWAXPY - Computes `w = alpha x + y`.
830: Logically Collective
832: Input Parameters:
833: + alpha - the scalar
834: . x - first vector, multiplied by `alpha`
835: - y - second vector
837: Output Parameter:
838: . w - the result
840: Level: intermediate
842: Note:
843: `w` cannot be either `x` or `y`, but `x` and `y` can be the same
845: Developer Notes:
846: The implementation is optimized for alpha of -1.0, 0.0, and 1.0
848: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
849: @*/
850: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
851: {
852: PetscFunctionBegin;
853: PetscCall(VecWAXPYAsync_Private(w, alpha, x, y, NULL));
854: PetscFunctionReturn(PETSC_SUCCESS);
855: }
857: /*@C
858: VecSetValues - Inserts or adds values into certain locations of a vector.
860: Not Collective
862: Input Parameters:
863: + x - vector to insert in
864: . ni - number of elements to add
865: . ix - indices where to add
866: . y - array of values
867: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries
869: Level: beginner
871: Notes:
872: .vb
873: `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
874: .ve
876: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
877: options cannot be mixed without intervening calls to the assembly
878: routines.
880: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
881: MUST be called after all calls to `VecSetValues()` have been completed.
883: VecSetValues() uses 0-based indices in Fortran as well as in C.
885: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
886: negative indices may be passed in ix. These rows are
887: simply ignored. This allows easily inserting element load matrices
888: with homogeneous Dirichlet boundary conditions that you don't want represented
889: in the vector.
891: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
892: `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
893: @*/
894: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
895: {
896: PetscFunctionBeginHot;
898: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
899: PetscAssertPointer(ix, 3);
900: PetscAssertPointer(y, 4);
903: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
904: PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
905: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
906: PetscCall(PetscObjectStateIncrease((PetscObject)x));
907: PetscFunctionReturn(PETSC_SUCCESS);
908: }
910: /*@C
911: VecGetValues - Gets values from certain locations of a vector. Currently
912: can only get values on the same processor on which they are owned
914: Not Collective
916: Input Parameters:
917: + x - vector to get values from
918: . ni - number of elements to get
919: - ix - indices where to get them from (in global 1d numbering)
921: Output Parameter:
922: . y - array of values
924: Level: beginner
926: Notes:
927: The user provides the allocated array y; it is NOT allocated in this routine
929: `VecGetValues()` gets y[i] = x[ix[i]], for i=0,...,ni-1.
931: `VecAssemblyBegin()` and `VecAssemblyEnd()` MUST be called before calling this if `VecSetValues()` or related routine has been called
933: VecGetValues() uses 0-based indices in Fortran as well as in C.
935: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
936: negative indices may be passed in ix. These rows are
937: simply ignored.
939: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
940: @*/
941: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
942: {
943: PetscFunctionBegin;
945: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
946: PetscAssertPointer(ix, 3);
947: PetscAssertPointer(y, 4);
949: VecCheckAssembled(x);
950: PetscUseTypeMethod(x, getvalues, ni, ix, y);
951: PetscFunctionReturn(PETSC_SUCCESS);
952: }
954: /*@C
955: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
957: Not Collective
959: Input Parameters:
960: + x - vector to insert in
961: . ni - number of blocks to add
962: . ix - indices where to add in block count, rather than element count
963: . y - array of values
964: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries
966: Level: intermediate
968: Notes:
969: `VecSetValuesBlocked()` sets x[bs*ix[i]+j] = y[bs*i+j],
970: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
972: Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
973: options cannot be mixed without intervening calls to the assembly
974: routines.
976: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
977: MUST be called after all calls to `VecSetValuesBlocked()` have been completed.
979: `VecSetValuesBlocked()` uses 0-based indices in Fortran as well as in C.
981: Negative indices may be passed in ix, these rows are
982: simply ignored. This allows easily inserting element load matrices
983: with homogeneous Dirichlet boundary conditions that you don't want represented
984: in the vector.
986: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
987: `VecSetValues()`
988: @*/
989: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
990: {
991: PetscFunctionBeginHot;
993: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
994: PetscAssertPointer(ix, 3);
995: PetscAssertPointer(y, 4);
998: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
999: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1000: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1001: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1002: PetscFunctionReturn(PETSC_SUCCESS);
1003: }
1005: /*@C
1006: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1007: using a local ordering of the nodes.
1009: Not Collective
1011: Input Parameters:
1012: + x - vector to insert in
1013: . ni - number of elements to add
1014: . ix - indices where to add
1015: . y - array of values
1016: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1018: Level: intermediate
1020: Notes:
1021: `VecSetValuesLocal()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
1023: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
1024: options cannot be mixed without intervening calls to the assembly
1025: routines.
1027: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1028: MUST be called after all calls to `VecSetValuesLocal()` have been completed.
1030: `VecSetValuesLocal()` uses 0-based indices in Fortran as well as in C.
1032: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1033: `VecSetValuesBlockedLocal()`
1034: @*/
1035: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1036: {
1037: PetscInt lixp[128], *lix = lixp;
1039: PetscFunctionBeginHot;
1041: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1042: PetscAssertPointer(ix, 3);
1043: PetscAssertPointer(y, 4);
1046: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1047: if (!x->ops->setvalueslocal) {
1048: if (x->map->mapping) {
1049: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1050: PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1051: PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1052: if (ni > 128) PetscCall(PetscFree(lix));
1053: } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1054: } else PetscUseTypeMethod(x, setvalueslocal, ni, ix, y, iora);
1055: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1056: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1057: PetscFunctionReturn(PETSC_SUCCESS);
1058: }
1060: /*@
1061: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1062: using a local ordering of the nodes.
1064: Not Collective
1066: Input Parameters:
1067: + x - vector to insert in
1068: . ni - number of blocks to add
1069: . ix - indices where to add in block count, not element count
1070: . y - array of values
1071: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1073: Level: intermediate
1075: Notes:
1076: `VecSetValuesBlockedLocal()` sets x[bs*ix[i]+j] = y[bs*i+j],
1077: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with `VecSetBlockSize()`.
1079: Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1080: options cannot be mixed without intervening calls to the assembly
1081: routines.
1083: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1084: MUST be called after all calls to `VecSetValuesBlockedLocal()` have been completed.
1086: `VecSetValuesBlockedLocal()` uses 0-based indices in Fortran as well as in C.
1088: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1089: `VecSetLocalToGlobalMapping()`
1090: @*/
1091: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1092: {
1093: PetscInt lixp[128], *lix = lixp;
1095: PetscFunctionBeginHot;
1097: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1098: PetscAssertPointer(ix, 3);
1099: PetscAssertPointer(y, 4);
1101: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1102: if (x->map->mapping) {
1103: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1104: PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1105: PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1106: if (ni > 128) PetscCall(PetscFree(lix));
1107: } else {
1108: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1109: }
1110: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1111: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1112: PetscFunctionReturn(PETSC_SUCCESS);
1113: }
1115: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1116: {
1117: PetscFunctionBegin;
1120: VecCheckAssembled(x);
1122: if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1123: PetscAssertPointer(y, 3);
1124: for (PetscInt i = 0; i < nv; ++i) {
1127: PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1128: VecCheckSameSize(x, 1, y[i], 3);
1129: VecCheckAssembled(y[i]);
1130: PetscCall(VecLockReadPush(y[i]));
1131: }
1132: PetscAssertPointer(result, 4);
1135: PetscCall(VecLockReadPush(x));
1136: PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1137: PetscCall((*mxdot)(x, nv, y, result));
1138: PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1139: PetscCall(VecLockReadPop(x));
1140: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1141: PetscFunctionReturn(PETSC_SUCCESS);
1142: }
1144: /*@
1145: VecMTDot - Computes indefinite vector multiple dot products.
1146: That is, it does NOT use the complex conjugate.
1148: Collective
1150: Input Parameters:
1151: + x - one vector
1152: . nv - number of vectors
1153: - y - array of vectors. Note that vectors are pointers
1155: Output Parameter:
1156: . val - array of the dot products
1158: Level: intermediate
1160: Notes for Users of Complex Numbers:
1161: For complex vectors, `VecMTDot()` computes the indefinite form
1162: $ val = (x,y) = y^T x,
1163: where y^T denotes the transpose of y.
1165: Use `VecMDot()` for the inner product
1166: $ val = (x,y) = y^H x,
1167: where y^H denotes the conjugate transpose of y.
1169: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1170: @*/
1171: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1172: {
1173: PetscFunctionBegin;
1175: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1176: PetscFunctionReturn(PETSC_SUCCESS);
1177: }
1179: /*@
1180: VecMDot - Computes multiple vector dot products.
1182: Collective
1184: Input Parameters:
1185: + x - one vector
1186: . nv - number of vectors
1187: - y - array of vectors.
1189: Output Parameter:
1190: . val - array of the dot products (does not allocate the array)
1192: Level: intermediate
1194: Notes for Users of Complex Numbers:
1195: For complex vectors, `VecMDot()` computes
1196: $ val = (x,y) = y^H x,
1197: where y^H denotes the conjugate transpose of y.
1199: Use `VecMTDot()` for the indefinite form
1200: $ val = (x,y) = y^T x,
1201: where y^T denotes the transpose of y.
1203: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`
1204: @*/
1205: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1206: {
1207: PetscFunctionBegin;
1209: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1210: PetscFunctionReturn(PETSC_SUCCESS);
1211: }
1213: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
1214: {
1215: PetscFunctionBegin;
1217: VecCheckAssembled(y);
1219: PetscCall(VecSetErrorIfLocked(y, 1));
1220: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1221: if (nv) {
1222: PetscInt zeros = 0;
1224: PetscAssertPointer(alpha, 3);
1225: PetscAssertPointer(x, 4);
1226: for (PetscInt i = 0; i < nv; ++i) {
1230: PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1231: VecCheckSameSize(y, 1, x[i], 4);
1232: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1233: VecCheckAssembled(x[i]);
1234: PetscCall(VecLockReadPush(x[i]));
1235: zeros += alpha[i] == (PetscScalar)0.0;
1236: }
1238: if (zeros < nv) {
1239: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1240: VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1241: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1242: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1243: }
1245: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1246: }
1247: PetscFunctionReturn(PETSC_SUCCESS);
1248: }
1250: /*@
1251: VecMAXPY - Computes `y = y + sum alpha[i] x[i]`
1253: Logically Collective
1255: Input Parameters:
1256: + nv - number of scalars and x-vectors
1257: . alpha - array of scalars
1258: . y - one vector
1259: - x - array of vectors
1261: Level: intermediate
1263: Note:
1264: `y` cannot be any of the `x` vectors
1266: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`,`VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1267: @*/
1268: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1269: {
1270: PetscFunctionBegin;
1271: PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1272: PetscFunctionReturn(PETSC_SUCCESS);
1273: }
1275: /*@
1276: VecMAXPBY - Computes `y = beta y + sum alpha[i] x[i]`
1278: Logically Collective
1280: Input Parameters:
1281: + nv - number of scalars and x-vectors
1282: . alpha - array of scalars
1283: . beta - scalar
1284: . y - one vector
1285: - x - array of vectors
1287: Level: intermediate
1289: Note:
1290: `y` cannot be any of the `x` vectors.
1292: Developer Notes:
1293: This is a convenience routine, but implementations might be able to optimize it, for example, when `beta` is zero.
1295: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1296: @*/
1297: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
1298: {
1299: PetscFunctionBegin;
1301: VecCheckAssembled(y);
1303: PetscCall(VecSetErrorIfLocked(y, 1));
1304: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1307: if (y->ops->maxpby) {
1308: PetscInt zeros = 0;
1310: if (nv) {
1311: PetscAssertPointer(alpha, 3);
1312: PetscAssertPointer(x, 5);
1313: }
1315: for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
1319: PetscCheckSameTypeAndComm(y, 1, x[i], 5);
1320: VecCheckSameSize(y, 1, x[i], 5);
1321: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1322: VecCheckAssembled(x[i]);
1323: PetscCall(VecLockReadPush(x[i]));
1324: zeros += alpha[i] == (PetscScalar)0.0;
1325: }
1327: if (zeros < nv) { // has nonzero alpha
1328: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1329: PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1330: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1331: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1332: } else {
1333: PetscCall(VecScale(y, beta));
1334: }
1336: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1337: } else { // no maxpby
1338: if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1339: else PetscCall(VecScale(y, beta));
1340: PetscCall(VecMAXPY(y, nv, alpha, x));
1341: }
1342: PetscFunctionReturn(PETSC_SUCCESS);
1343: }
1345: /*@
1346: VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
1347: in the order they appear in the array. The concatenated vector resides on the same
1348: communicator and is the same type as the source vectors.
1350: Collective
1352: Input Parameters:
1353: + nx - number of vectors to be concatenated
1354: - X - array containing the vectors to be concatenated in the order of concatenation
1356: Output Parameters:
1357: + Y - concatenated vector
1358: - x_is - array of index sets corresponding to the concatenated components of `Y` (pass `NULL` if not needed)
1360: Level: advanced
1362: Notes:
1363: Concatenation is similar to the functionality of a `VECNEST` object; they both represent combination of
1364: different vector spaces. However, concatenated vectors do not store any information about their
1365: sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1366: manipulation of data in the concatenated vector that corresponds to the original components at creation.
1368: This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1369: has to operate on combined vector spaces and cannot utilize `VECNEST` objects due to incompatibility with
1370: bound projections.
1372: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1373: @*/
1374: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1375: {
1376: MPI_Comm comm;
1377: VecType vec_type;
1378: Vec Ytmp, Xtmp;
1379: IS *is_tmp;
1380: PetscInt i, shift = 0, Xnl, Xng, Xbegin;
1382: PetscFunctionBegin;
1386: PetscAssertPointer(Y, 3);
1388: if ((*X)->ops->concatenate) {
1389: /* use the dedicated concatenation function if available */
1390: PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1391: } else {
1392: /* loop over vectors and start creating IS */
1393: comm = PetscObjectComm((PetscObject)(*X));
1394: PetscCall(VecGetType(*X, &vec_type));
1395: PetscCall(PetscMalloc1(nx, &is_tmp));
1396: for (i = 0; i < nx; i++) {
1397: PetscCall(VecGetSize(X[i], &Xng));
1398: PetscCall(VecGetLocalSize(X[i], &Xnl));
1399: PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1400: PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1401: shift += Xng;
1402: }
1403: /* create the concatenated vector */
1404: PetscCall(VecCreate(comm, &Ytmp));
1405: PetscCall(VecSetType(Ytmp, vec_type));
1406: PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1407: PetscCall(VecSetUp(Ytmp));
1408: /* copy data from X array to Y and return */
1409: for (i = 0; i < nx; i++) {
1410: PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1411: PetscCall(VecCopy(X[i], Xtmp));
1412: PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1413: }
1414: *Y = Ytmp;
1415: if (x_is) {
1416: *x_is = is_tmp;
1417: } else {
1418: for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1419: PetscCall(PetscFree(is_tmp));
1420: }
1421: }
1422: PetscFunctionReturn(PETSC_SUCCESS);
1423: }
1425: /* A helper function for VecGetSubVector to check if we can implement it with no-copy (i.e. the subvector shares
1426: memory with the original vector), and the block size of the subvector.
1428: Input Parameters:
1429: + X - the original vector
1430: - is - the index set of the subvector
1432: Output Parameters:
1433: + contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
1434: . start - start of contiguous block, as an offset from the start of the ownership range of the original vector
1435: - blocksize - the block size of the subvector
1437: */
1438: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1439: {
1440: PetscInt gstart, gend, lstart;
1441: PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1442: PetscInt n, N, ibs, vbs, bs = -1;
1444: PetscFunctionBegin;
1445: PetscCall(ISGetLocalSize(is, &n));
1446: PetscCall(ISGetSize(is, &N));
1447: PetscCall(ISGetBlockSize(is, &ibs));
1448: PetscCall(VecGetBlockSize(X, &vbs));
1449: PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1450: PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1451: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1452: if (ibs > 1) {
1453: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1454: bs = ibs;
1455: } else {
1456: if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1457: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1458: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1459: }
1461: *contig = red[0];
1462: *start = lstart;
1463: *blocksize = bs;
1464: PetscFunctionReturn(PETSC_SUCCESS);
1465: }
1467: /* A helper function for VecGetSubVector, to be used when we have to build a standalone subvector through VecScatter
1469: Input Parameters:
1470: + X - the original vector
1471: . is - the index set of the subvector
1472: - bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()
1474: Output Parameter:
1475: . Z - the subvector, which will compose the VecScatter context on output
1476: */
1477: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1478: {
1479: PetscInt n, N;
1480: VecScatter vscat;
1481: Vec Y;
1483: PetscFunctionBegin;
1484: PetscCall(ISGetLocalSize(is, &n));
1485: PetscCall(ISGetSize(is, &N));
1486: PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1487: PetscCall(VecSetSizes(Y, n, N));
1488: PetscCall(VecSetBlockSize(Y, bs));
1489: PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1490: PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1491: PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1492: PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1493: PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1494: PetscCall(VecScatterDestroy(&vscat));
1495: *Z = Y;
1496: PetscFunctionReturn(PETSC_SUCCESS);
1497: }
1499: /*@
1500: VecGetSubVector - Gets a vector representing part of another vector
1502: Collective
1504: Input Parameters:
1505: + X - vector from which to extract a subvector
1506: - is - index set representing portion of `X` to extract
1508: Output Parameter:
1509: . Y - subvector corresponding to `is`
1511: Level: advanced
1513: Notes:
1514: The subvector `Y` should be returned with `VecRestoreSubVector()`.
1515: `X` and must be defined on the same communicator
1517: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1518: modifying the subvector. Other non-overlapping subvectors can still be obtained from X using this function.
1520: 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`.
1522: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1523: @*/
1524: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1525: {
1526: Vec Z;
1528: PetscFunctionBegin;
1531: PetscCheckSameComm(X, 1, is, 2);
1532: PetscAssertPointer(Y, 3);
1533: if (X->ops->getsubvector) {
1534: PetscUseTypeMethod(X, getsubvector, is, &Z);
1535: } else { /* Default implementation currently does no caching */
1536: PetscBool contig;
1537: PetscInt n, N, start, bs;
1539: PetscCall(ISGetLocalSize(is, &n));
1540: PetscCall(ISGetSize(is, &N));
1541: PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1542: if (contig) { /* We can do a no-copy implementation */
1543: const PetscScalar *x;
1544: PetscInt state = 0;
1545: PetscBool isstd, iscuda, iship;
1547: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1548: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1549: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1550: if (iscuda) {
1551: #if defined(PETSC_HAVE_CUDA)
1552: const PetscScalar *x_d;
1553: PetscMPIInt size;
1554: PetscOffloadMask flg;
1556: PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1557: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1558: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1559: if (x) x += start;
1560: if (x_d) x_d += start;
1561: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1562: if (size == 1) {
1563: PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1564: } else {
1565: PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1566: }
1567: Z->offloadmask = flg;
1568: #endif
1569: } else if (iship) {
1570: #if defined(PETSC_HAVE_HIP)
1571: const PetscScalar *x_d;
1572: PetscMPIInt size;
1573: PetscOffloadMask flg;
1575: PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1576: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1577: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1578: if (x) x += start;
1579: if (x_d) x_d += start;
1580: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1581: if (size == 1) {
1582: PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1583: } else {
1584: PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1585: }
1586: Z->offloadmask = flg;
1587: #endif
1588: } else if (isstd) {
1589: PetscMPIInt size;
1591: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1592: PetscCall(VecGetArrayRead(X, &x));
1593: if (x) x += start;
1594: if (size == 1) {
1595: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1596: } else {
1597: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1598: }
1599: PetscCall(VecRestoreArrayRead(X, &x));
1600: } else { /* default implementation: use place array */
1601: PetscCall(VecGetArrayRead(X, &x));
1602: PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1603: PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1604: PetscCall(VecSetSizes(Z, n, N));
1605: PetscCall(VecSetBlockSize(Z, bs));
1606: PetscCall(VecPlaceArray(Z, x ? x + start : NULL));
1607: PetscCall(VecRestoreArrayRead(X, &x));
1608: }
1610: /* this is relevant only in debug mode */
1611: PetscCall(VecLockGet(X, &state));
1612: if (state) PetscCall(VecLockReadPush(Z));
1613: Z->ops->placearray = NULL;
1614: Z->ops->replacearray = NULL;
1615: } else { /* Have to create a scatter and do a copy */
1616: PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1617: }
1618: }
1619: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1620: if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1621: PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1622: *Y = Z;
1623: PetscFunctionReturn(PETSC_SUCCESS);
1624: }
1626: /*@
1627: VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`
1629: Collective
1631: Input Parameters:
1632: + X - vector from which subvector was obtained
1633: . is - index set representing the subset of `X`
1634: - Y - subvector being restored
1636: Level: advanced
1638: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1639: @*/
1640: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1641: {
1642: PETSC_UNUSED PetscObjectState dummystate = 0;
1643: PetscBool unchanged;
1645: PetscFunctionBegin;
1648: PetscCheckSameComm(X, 1, is, 2);
1649: PetscAssertPointer(Y, 3);
1652: if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1653: else {
1654: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1655: if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1656: VecScatter scatter;
1657: PetscInt state;
1659: PetscCall(VecLockGet(X, &state));
1660: PetscCheck(state == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec X is locked for read-only or read/write access");
1662: PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1663: if (scatter) {
1664: PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1665: PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1666: } else {
1667: PetscBool iscuda, iship;
1668: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1669: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1671: if (iscuda) {
1672: #if defined(PETSC_HAVE_CUDA)
1673: PetscOffloadMask ymask = (*Y)->offloadmask;
1675: /* The offloadmask of X dictates where to move memory
1676: If X GPU data is valid, then move Y data on GPU if needed
1677: Otherwise, move back to the CPU */
1678: switch (X->offloadmask) {
1679: case PETSC_OFFLOAD_BOTH:
1680: if (ymask == PETSC_OFFLOAD_CPU) {
1681: PetscCall(VecCUDAResetArray(*Y));
1682: } else if (ymask == PETSC_OFFLOAD_GPU) {
1683: X->offloadmask = PETSC_OFFLOAD_GPU;
1684: }
1685: break;
1686: case PETSC_OFFLOAD_GPU:
1687: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1688: break;
1689: case PETSC_OFFLOAD_CPU:
1690: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1691: break;
1692: case PETSC_OFFLOAD_UNALLOCATED:
1693: case PETSC_OFFLOAD_KOKKOS:
1694: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1695: }
1696: #endif
1697: } else if (iship) {
1698: #if defined(PETSC_HAVE_HIP)
1699: PetscOffloadMask ymask = (*Y)->offloadmask;
1701: /* The offloadmask of X dictates where to move memory
1702: If X GPU data is valid, then move Y data on GPU if needed
1703: Otherwise, move back to the CPU */
1704: switch (X->offloadmask) {
1705: case PETSC_OFFLOAD_BOTH:
1706: if (ymask == PETSC_OFFLOAD_CPU) {
1707: PetscCall(VecHIPResetArray(*Y));
1708: } else if (ymask == PETSC_OFFLOAD_GPU) {
1709: X->offloadmask = PETSC_OFFLOAD_GPU;
1710: }
1711: break;
1712: case PETSC_OFFLOAD_GPU:
1713: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1714: break;
1715: case PETSC_OFFLOAD_CPU:
1716: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1717: break;
1718: case PETSC_OFFLOAD_UNALLOCATED:
1719: case PETSC_OFFLOAD_KOKKOS:
1720: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1721: }
1722: #endif
1723: } else {
1724: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1725: PetscCall(VecResetArray(*Y));
1726: }
1727: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1728: }
1729: }
1730: }
1731: PetscCall(VecDestroy(Y));
1732: PetscFunctionReturn(PETSC_SUCCESS);
1733: }
1735: /*@
1736: VecCreateLocalVector - Creates a vector object suitable for use with `VecGetLocalVector()` and friends. You must call `VecDestroy()` when the
1737: vector is no longer needed.
1739: Not Collective.
1741: Input Parameter:
1742: . v - The vector for which the local vector is desired.
1744: Output Parameter:
1745: . w - Upon exit this contains the local vector.
1747: Level: beginner
1749: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1750: @*/
1751: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1752: {
1753: PetscMPIInt size;
1755: PetscFunctionBegin;
1757: PetscAssertPointer(w, 2);
1758: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1759: if (size == 1) PetscCall(VecDuplicate(v, w));
1760: else if (v->ops->createlocalvector) PetscUseTypeMethod(v, createlocalvector, w);
1761: else {
1762: VecType type;
1763: PetscInt n;
1765: PetscCall(VecCreate(PETSC_COMM_SELF, w));
1766: PetscCall(VecGetLocalSize(v, &n));
1767: PetscCall(VecSetSizes(*w, n, n));
1768: PetscCall(VecGetBlockSize(v, &n));
1769: PetscCall(VecSetBlockSize(*w, n));
1770: PetscCall(VecGetType(v, &type));
1771: PetscCall(VecSetType(*w, type));
1772: }
1773: PetscFunctionReturn(PETSC_SUCCESS);
1774: }
1776: /*@
1777: VecGetLocalVectorRead - Maps the local portion of a vector into a
1778: vector.
1780: Not Collective.
1782: Input Parameter:
1783: . v - The vector for which the local vector is desired.
1785: Output Parameter:
1786: . w - Upon exit this contains the local vector.
1788: Level: beginner
1790: Notes:
1791: You must call `VecRestoreLocalVectorRead()` when the local
1792: vector is no longer needed.
1794: This function is similar to `VecGetArrayRead()` which maps the local
1795: portion into a raw pointer. `VecGetLocalVectorRead()` is usually
1796: almost as efficient as `VecGetArrayRead()` but in certain circumstances
1797: `VecGetLocalVectorRead()` can be much more efficient than
1798: `VecGetArrayRead()`. This is because the construction of a contiguous
1799: array representing the vector data required by `VecGetArrayRead()` can
1800: be an expensive operation for certain vector types. For example, for
1801: GPU vectors `VecGetArrayRead()` requires that the data between device
1802: and host is synchronized.
1804: Unlike `VecGetLocalVector()`, this routine is not collective and
1805: preserves cached information.
1807: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1808: @*/
1809: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1810: {
1811: PetscFunctionBegin;
1814: VecCheckSameLocalSize(v, 1, w, 2);
1815: if (v->ops->getlocalvectorread) {
1816: PetscUseTypeMethod(v, getlocalvectorread, w);
1817: } else {
1818: PetscScalar *a;
1820: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1821: PetscCall(VecPlaceArray(w, a));
1822: }
1823: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1824: PetscCall(VecLockReadPush(v));
1825: PetscCall(VecLockReadPush(w));
1826: PetscFunctionReturn(PETSC_SUCCESS);
1827: }
1829: /*@
1830: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1831: previously mapped into a vector using `VecGetLocalVectorRead()`.
1833: Not Collective.
1835: Input Parameters:
1836: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
1837: - w - The vector into which the local portion of `v` was mapped.
1839: Level: beginner
1841: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1842: @*/
1843: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1844: {
1845: PetscFunctionBegin;
1848: if (v->ops->restorelocalvectorread) {
1849: PetscUseTypeMethod(v, restorelocalvectorread, w);
1850: } else {
1851: const PetscScalar *a;
1853: PetscCall(VecGetArrayRead(w, &a));
1854: PetscCall(VecRestoreArrayRead(v, &a));
1855: PetscCall(VecResetArray(w));
1856: }
1857: PetscCall(VecLockReadPop(v));
1858: PetscCall(VecLockReadPop(w));
1859: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1860: PetscFunctionReturn(PETSC_SUCCESS);
1861: }
1863: /*@
1864: VecGetLocalVector - Maps the local portion of a vector into a
1865: vector.
1867: Collective
1869: Input Parameter:
1870: . v - The vector for which the local vector is desired.
1872: Output Parameter:
1873: . w - Upon exit this contains the local vector.
1875: Level: beginner
1877: Notes:
1878: You must call `VecRestoreLocalVector()` when the local
1879: vector is no longer needed.
1881: This function is similar to `VecGetArray()` which maps the local
1882: portion into a raw pointer. `VecGetLocalVector()` is usually about as
1883: efficient as `VecGetArray()` but in certain circumstances
1884: `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1885: This is because the construction of a contiguous array representing
1886: the vector data required by `VecGetArray()` can be an expensive
1887: operation for certain vector types. For example, for GPU vectors
1888: `VecGetArray()` requires that the data between device and host is
1889: synchronized.
1891: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1892: @*/
1893: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1894: {
1895: PetscFunctionBegin;
1898: VecCheckSameLocalSize(v, 1, w, 2);
1899: if (v->ops->getlocalvector) {
1900: PetscUseTypeMethod(v, getlocalvector, w);
1901: } else {
1902: PetscScalar *a;
1904: PetscCall(VecGetArray(v, &a));
1905: PetscCall(VecPlaceArray(w, a));
1906: }
1907: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1908: PetscFunctionReturn(PETSC_SUCCESS);
1909: }
1911: /*@
1912: VecRestoreLocalVector - Unmaps the local portion of a vector
1913: previously mapped into a vector using `VecGetLocalVector()`.
1915: Logically Collective.
1917: Input Parameters:
1918: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
1919: - w - The vector into which the local portion of `v` was mapped.
1921: Level: beginner
1923: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1924: @*/
1925: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1926: {
1927: PetscFunctionBegin;
1930: if (v->ops->restorelocalvector) {
1931: PetscUseTypeMethod(v, restorelocalvector, w);
1932: } else {
1933: PetscScalar *a;
1934: PetscCall(VecGetArray(w, &a));
1935: PetscCall(VecRestoreArray(v, &a));
1936: PetscCall(VecResetArray(w));
1937: }
1938: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1939: PetscCall(PetscObjectStateIncrease((PetscObject)v));
1940: PetscFunctionReturn(PETSC_SUCCESS);
1941: }
1943: /*@C
1944: VecGetArray - Returns a pointer to a contiguous array that contains this
1945: MPI processes's portion of the vector data
1947: Logically Collective
1949: Input Parameter:
1950: . x - the vector
1952: Output Parameter:
1953: . a - location to put pointer to the array
1955: Level: beginner
1957: Notes:
1958: For the standard PETSc vectors, `VecGetArray()` returns a pointer to the local data array and
1959: does not use any copies. If the underlying vector data is not stored in a contiguous array
1960: this routine will copy the data to a contiguous array and return a pointer to that. You MUST
1961: call `VecRestoreArray()` when you no longer need access to the array.
1963: Fortran Notes:
1964: `VecGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayF90()`
1966: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
1967: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
1968: @*/
1969: PetscErrorCode VecGetArray(Vec x, PetscScalar **a)
1970: {
1971: PetscFunctionBegin;
1973: PetscCall(VecSetErrorIfLocked(x, 1));
1974: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
1975: PetscUseTypeMethod(x, getarray, a);
1976: } else if (x->petscnative) { /* VECSTANDARD */
1977: *a = *((PetscScalar **)x->data);
1978: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
1979: PetscFunctionReturn(PETSC_SUCCESS);
1980: }
1982: /*@C
1983: VecRestoreArray - Restores a vector after `VecGetArray()` has been called and the array is no longer needed
1985: Logically Collective
1987: Input Parameters:
1988: + x - the vector
1989: - a - location of pointer to array obtained from `VecGetArray()`
1991: Level: beginner
1993: Fortran Notes:
1994: `VecRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayF90()`
1996: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
1997: `VecGetArrayPair()`, `VecRestoreArrayPair()`
1998: @*/
1999: PetscErrorCode VecRestoreArray(Vec x, PetscScalar **a)
2000: {
2001: PetscFunctionBegin;
2003: if (a) PetscAssertPointer(a, 2);
2004: if (x->ops->restorearray) {
2005: PetscUseTypeMethod(x, restorearray, a);
2006: } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
2007: if (a) *a = NULL;
2008: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2009: PetscFunctionReturn(PETSC_SUCCESS);
2010: }
2011: /*@C
2012: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
2014: Not Collective
2016: Input Parameter:
2017: . x - the vector
2019: Output Parameter:
2020: . a - the array
2022: Level: beginner
2024: Notes:
2025: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2027: Unlike `VecGetArray()`, preserves cached information like vector norms.
2029: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
2030: implementations may require a copy, but such implementations should cache the contiguous representation so that
2031: only one copy is performed when this routine is called multiple times in sequence.
2033: Fortran Notes:
2034: `VecGetArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayReadF90()`
2036: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2037: @*/
2038: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar **a)
2039: {
2040: PetscFunctionBegin;
2042: PetscAssertPointer(a, 2);
2043: if (x->ops->getarrayread) {
2044: PetscUseTypeMethod(x, getarrayread, a);
2045: } else if (x->ops->getarray) {
2046: PetscObjectState state;
2048: /* VECNEST, VECCUDA, VECKOKKOS etc */
2049: // x->ops->getarray may bump the object state, but since we know this is a read-only get
2050: // we can just undo that
2051: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2052: PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
2053: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2054: } else if (x->petscnative) {
2055: /* VECSTANDARD */
2056: *a = *((PetscScalar **)x->data);
2057: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
2058: PetscFunctionReturn(PETSC_SUCCESS);
2059: }
2061: /*@C
2062: VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`
2064: Not Collective
2066: Input Parameters:
2067: + x - the vector
2068: - a - the array
2070: Level: beginner
2072: Fortran Notes:
2073: `VecRestoreArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayReadF90()`
2075: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2076: @*/
2077: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar **a)
2078: {
2079: PetscFunctionBegin;
2081: if (a) PetscAssertPointer(a, 2);
2082: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
2083: /* nothing */
2084: } else if (x->ops->restorearrayread) { /* VECNEST */
2085: PetscUseTypeMethod(x, restorearrayread, a);
2086: } else { /* No one? */
2087: PetscObjectState state;
2089: // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2090: // we can just undo that
2091: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2092: PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2093: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2094: }
2095: if (a) *a = NULL;
2096: PetscFunctionReturn(PETSC_SUCCESS);
2097: }
2099: /*@C
2100: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
2101: MPI processes's portion of the vector data.
2103: Logically Collective
2105: Input Parameter:
2106: . x - the vector
2108: Output Parameter:
2109: . a - location to put pointer to the array
2111: Level: intermediate
2113: Note:
2114: The values in this array are NOT valid, the caller of this routine is responsible for putting
2115: values into the array; any values it does not set will be invalid.
2117: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2119: For vectors associated with GPUs, the host and device vectors are not synchronized before
2120: giving access. If you need correct values in the array use `VecGetArray()`
2122: Fortran Notes:
2123: `VecGetArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayWriteF90()`
2125: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteF90()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
2126: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
2127: @*/
2128: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar **a)
2129: {
2130: PetscFunctionBegin;
2132: PetscAssertPointer(a, 2);
2133: PetscCall(VecSetErrorIfLocked(x, 1));
2134: if (x->ops->getarraywrite) {
2135: PetscUseTypeMethod(x, getarraywrite, a);
2136: } else {
2137: PetscCall(VecGetArray(x, a));
2138: }
2139: PetscFunctionReturn(PETSC_SUCCESS);
2140: }
2142: /*@C
2143: VecRestoreArrayWrite - Restores a vector after `VecGetArrayWrite()` has been called.
2145: Logically Collective
2147: Input Parameters:
2148: + x - the vector
2149: - a - location of pointer to array obtained from `VecGetArray()`
2151: Level: beginner
2153: Fortran Notes:
2154: `VecRestoreArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayWriteF90()`
2156: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteF90()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2157: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2158: @*/
2159: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar **a)
2160: {
2161: PetscFunctionBegin;
2163: if (a) PetscAssertPointer(a, 2);
2164: if (x->ops->restorearraywrite) {
2165: PetscUseTypeMethod(x, restorearraywrite, a);
2166: } else if (x->ops->restorearray) {
2167: PetscUseTypeMethod(x, restorearray, a);
2168: }
2169: if (a) *a = NULL;
2170: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2171: PetscFunctionReturn(PETSC_SUCCESS);
2172: }
2174: /*@C
2175: VecGetArrays - Returns a pointer to the arrays in a set of vectors
2176: that were created by a call to `VecDuplicateVecs()`.
2178: Logically Collective; No Fortran Support
2180: Input Parameters:
2181: + x - the vectors
2182: - n - the number of vectors
2184: Output Parameter:
2185: . a - location to put pointer to the array
2187: Level: intermediate
2189: Note:
2190: You MUST call `VecRestoreArrays()` when you no longer need access to the arrays.
2192: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2193: @*/
2194: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2195: {
2196: PetscInt i;
2197: PetscScalar **q;
2199: PetscFunctionBegin;
2200: PetscAssertPointer(x, 1);
2202: PetscAssertPointer(a, 3);
2203: PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2204: PetscCall(PetscMalloc1(n, &q));
2205: for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2206: *a = q;
2207: PetscFunctionReturn(PETSC_SUCCESS);
2208: }
2210: /*@C
2211: VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2212: has been called.
2214: Logically Collective; No Fortran Support
2216: Input Parameters:
2217: + x - the vector
2218: . n - the number of vectors
2219: - a - location of pointer to arrays obtained from `VecGetArrays()`
2221: Notes:
2222: For regular PETSc vectors this routine does not involve any copies. For
2223: any special vectors that do not store local vector data in a contiguous
2224: array, this routine will copy the data back into the underlying
2225: vector data structure from the arrays obtained with `VecGetArrays()`.
2227: Level: intermediate
2229: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2230: @*/
2231: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2232: {
2233: PetscInt i;
2234: PetscScalar **q = *a;
2236: PetscFunctionBegin;
2237: PetscAssertPointer(x, 1);
2239: PetscAssertPointer(a, 3);
2241: for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2242: PetscCall(PetscFree(q));
2243: PetscFunctionReturn(PETSC_SUCCESS);
2244: }
2246: /*@C
2247: VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g.,
2248: `VECCUDA`), the returned pointer will be a device pointer to the device memory that contains
2249: this MPI processes's portion of the vector data.
2251: Logically Collective; No Fortran Support
2253: Input Parameter:
2254: . x - the vector
2256: Output Parameters:
2257: + a - location to put pointer to the array
2258: - mtype - memory type of the array
2260: Level: beginner
2262: Note:
2263: Device data is guaranteed to have the latest value. Otherwise, when this is a host vector
2264: (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host
2265: pointer.
2267: For `VECKOKKOS`, if Kokkos is configured without device (e.g., use serial or openmp), per
2268: this function, the vector works like `VECSEQ`/`VECMPI`; otherwise, it works like `VECCUDA` or
2269: `VECHIP` etc.
2271: Use `VecRestoreArrayAndMemType()` when the array access is no longer needed.
2273: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`,
2274: `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2275: @*/
2276: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2277: {
2278: PetscFunctionBegin;
2281: PetscAssertPointer(a, 2);
2282: if (mtype) PetscAssertPointer(mtype, 3);
2283: PetscCall(VecSetErrorIfLocked(x, 1));
2284: if (x->ops->getarrayandmemtype) {
2285: /* VECCUDA, VECKOKKOS etc */
2286: PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2287: } else {
2288: /* VECSTANDARD, VECNEST, VECVIENNACL */
2289: PetscCall(VecGetArray(x, a));
2290: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2291: }
2292: PetscFunctionReturn(PETSC_SUCCESS);
2293: }
2295: /*@C
2296: VecRestoreArrayAndMemType - Restores a vector after `VecGetArrayAndMemType()` has been called.
2298: Logically Collective; No Fortran Support
2300: Input Parameters:
2301: + x - the vector
2302: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`
2304: Level: beginner
2306: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`,
2307: `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2308: @*/
2309: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar **a)
2310: {
2311: PetscFunctionBegin;
2314: if (a) PetscAssertPointer(a, 2);
2315: if (x->ops->restorearrayandmemtype) {
2316: /* VECCUDA, VECKOKKOS etc */
2317: PetscUseTypeMethod(x, restorearrayandmemtype, a);
2318: } else {
2319: /* VECNEST, VECVIENNACL */
2320: PetscCall(VecRestoreArray(x, a));
2321: } /* VECSTANDARD does nothing */
2322: if (a) *a = NULL;
2323: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2324: PetscFunctionReturn(PETSC_SUCCESS);
2325: }
2327: /*@C
2328: VecGetArrayReadAndMemType - Like `VecGetArrayRead()`, but if the input vector is a device vector, it will return a read-only device pointer.
2329: The returned pointer is guaranteed to point to up-to-date data. For host vectors, it functions as `VecGetArrayRead()`.
2331: Not Collective; No Fortran Support
2333: Input Parameter:
2334: . x - the vector
2336: Output Parameters:
2337: + a - the array
2338: - mtype - memory type of the array
2340: Level: beginner
2342: Notes:
2343: The array must be returned using a matching call to `VecRestoreArrayReadAndMemType()`.
2345: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2346: @*/
2347: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar **a, PetscMemType *mtype)
2348: {
2349: PetscFunctionBegin;
2352: PetscAssertPointer(a, 2);
2353: if (mtype) PetscAssertPointer(mtype, 3);
2354: if (x->ops->getarrayreadandmemtype) {
2355: /* VECCUDA/VECHIP though they are also petscnative */
2356: PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2357: } else if (x->ops->getarrayandmemtype) {
2358: /* VECKOKKOS */
2359: PetscObjectState state;
2361: // see VecGetArrayRead() for why
2362: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2363: PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2364: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2365: } else {
2366: PetscCall(VecGetArrayRead(x, a));
2367: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2368: }
2369: PetscFunctionReturn(PETSC_SUCCESS);
2370: }
2372: /*@C
2373: VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`
2375: Not Collective; No Fortran Support
2377: Input Parameters:
2378: + x - the vector
2379: - a - the array
2381: Level: beginner
2383: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2384: @*/
2385: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar **a)
2386: {
2387: PetscFunctionBegin;
2390: if (a) PetscAssertPointer(a, 2);
2391: if (x->ops->restorearrayreadandmemtype) {
2392: /* VECCUDA/VECHIP */
2393: PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2394: } else if (!x->petscnative) {
2395: /* VECNEST */
2396: PetscCall(VecRestoreArrayRead(x, a));
2397: }
2398: if (a) *a = NULL;
2399: PetscFunctionReturn(PETSC_SUCCESS);
2400: }
2402: /*@C
2403: VecGetArrayWriteAndMemType - Like `VecGetArrayWrite()`, but if this is a device vector it will always return
2404: a device pointer to the device memory that contains this processor's portion of the vector data.
2406: Not Collective; No Fortran Support
2408: Input Parameter:
2409: . x - the vector
2411: Output Parameters:
2412: + a - the array
2413: - mtype - memory type of the array
2415: Level: beginner
2417: Note:
2418: The array must be returned using a matching call to `VecRestoreArrayWriteAndMemType()`, where it will label the device memory as most recent.
2420: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2421: @*/
2422: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2423: {
2424: PetscFunctionBegin;
2427: PetscCall(VecSetErrorIfLocked(x, 1));
2428: PetscAssertPointer(a, 2);
2429: if (mtype) PetscAssertPointer(mtype, 3);
2430: if (x->ops->getarraywriteandmemtype) {
2431: /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2432: PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2433: } else if (x->ops->getarrayandmemtype) {
2434: PetscCall(VecGetArrayAndMemType(x, a, mtype));
2435: } else {
2436: /* VECNEST, VECVIENNACL */
2437: PetscCall(VecGetArrayWrite(x, a));
2438: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2439: }
2440: PetscFunctionReturn(PETSC_SUCCESS);
2441: }
2443: /*@C
2444: VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`
2446: Not Collective; No Fortran Support
2448: Input Parameters:
2449: + x - the vector
2450: - a - the array
2452: Level: beginner
2454: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2455: @*/
2456: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar **a)
2457: {
2458: PetscFunctionBegin;
2461: PetscCall(VecSetErrorIfLocked(x, 1));
2462: if (a) PetscAssertPointer(a, 2);
2463: if (x->ops->restorearraywriteandmemtype) {
2464: /* VECCUDA/VECHIP */
2465: PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2466: PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2467: } else if (x->ops->restorearrayandmemtype) {
2468: PetscCall(VecRestoreArrayAndMemType(x, a));
2469: } else {
2470: PetscCall(VecRestoreArray(x, a));
2471: }
2472: if (a) *a = NULL;
2473: PetscFunctionReturn(PETSC_SUCCESS);
2474: }
2476: /*@
2477: VecPlaceArray - Allows one to replace the array in a vector with an
2478: array provided by the user. This is useful to avoid copying an array
2479: into a vector.
2481: Not Collective; No Fortran Support
2483: Input Parameters:
2484: + vec - the vector
2485: - array - the array
2487: Level: developer
2489: Notes:
2490: Use `VecReplaceArray()` instead to permanently replace the array
2492: You can return to the original array with a call to `VecResetArray()`. `vec` does not take
2493: ownership of `array` in any way.
2495: The user must free `array` themselves but be careful not to
2496: do so before the vector has either been destroyed, had its original array restored with
2497: `VecResetArray()` or permanently replaced with `VecReplaceArray()`.
2499: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2500: @*/
2501: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2502: {
2503: PetscFunctionBegin;
2506: if (array) PetscAssertPointer(array, 2);
2507: PetscUseTypeMethod(vec, placearray, array);
2508: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2509: PetscFunctionReturn(PETSC_SUCCESS);
2510: }
2512: /*@C
2513: VecReplaceArray - Allows one to replace the array in a vector with an
2514: array provided by the user. This is useful to avoid copying an array
2515: into a vector.
2517: Not Collective; No Fortran Support
2519: Input Parameters:
2520: + vec - the vector
2521: - array - the array
2523: Level: developer
2525: Notes:
2526: This permanently replaces the array and frees the memory associated
2527: with the old array. Use `VecPlaceArray()` to temporarily replace the array.
2529: The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
2530: freed by the user. It will be freed when the vector is destroyed.
2532: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2533: @*/
2534: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2535: {
2536: PetscFunctionBegin;
2539: PetscUseTypeMethod(vec, replacearray, array);
2540: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2541: PetscFunctionReturn(PETSC_SUCCESS);
2542: }
2544: /*MC
2545: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
2546: and makes them accessible via a Fortran pointer.
2548: Synopsis:
2549: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
2551: Collective
2553: Input Parameters:
2554: + x - a vector to mimic
2555: - n - the number of vectors to obtain
2557: Output Parameters:
2558: + y - Fortran pointer to the array of vectors
2559: - ierr - error code
2561: Example of Usage:
2562: .vb
2563: #include <petsc/finclude/petscvec.h>
2564: use petscvec
2566: Vec x
2567: Vec, pointer :: y(:)
2568: ....
2569: call VecDuplicateVecsF90(x,2,y,ierr)
2570: call VecSet(y(2),alpha,ierr)
2571: call VecSet(y(2),alpha,ierr)
2572: ....
2573: call VecDestroyVecsF90(2,y,ierr)
2574: .ve
2576: Level: beginner
2578: Note:
2579: Use `VecDestroyVecsF90()` to free the space.
2581: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecsF90()`, `VecDuplicateVecs()`
2582: M*/
2584: /*MC
2585: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2586: `VecGetArrayF90()`.
2588: Synopsis:
2589: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2591: Logically Collective
2593: Input Parameters:
2594: + x - vector
2595: - xx_v - the Fortran pointer to the array
2597: Output Parameter:
2598: . ierr - error code
2600: Example of Usage:
2601: .vb
2602: #include <petsc/finclude/petscvec.h>
2603: use petscvec
2605: PetscScalar, pointer :: xx_v(:)
2606: ....
2607: call VecGetArrayF90(x,xx_v,ierr)
2608: xx_v(3) = a
2609: call VecRestoreArrayF90(x,xx_v,ierr)
2610: .ve
2612: Level: beginner
2614: .seealso: [](ch_vectors), `Vec`, `VecGetArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrayReadF90()`
2615: M*/
2617: /*MC
2618: VecDestroyVecsF90 - Frees a block of vectors obtained with `VecDuplicateVecsF90()`.
2620: Synopsis:
2621: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
2623: Collective
2625: Input Parameters:
2626: + n - the number of vectors previously obtained
2627: - x - pointer to array of vector pointers
2629: Output Parameter:
2630: . ierr - error code
2632: Level: beginner
2634: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecs()`, `VecDuplicateVecsF90()`
2635: M*/
2637: /*MC
2638: VecGetArrayF90 - Accesses a vector array from Fortran. For default PETSc
2639: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2640: this routine is implementation dependent. You MUST call `VecRestoreArrayF90()`
2641: when you no longer need access to the array.
2643: Synopsis:
2644: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2646: Logically Collective
2648: Input Parameter:
2649: . x - vector
2651: Output Parameters:
2652: + xx_v - the Fortran pointer to the array
2653: - ierr - error code
2655: Example of Usage:
2656: .vb
2657: #include <petsc/finclude/petscvec.h>
2658: use petscvec
2660: PetscScalar, pointer :: xx_v(:)
2661: ....
2662: call VecGetArrayF90(x,xx_v,ierr)
2663: xx_v(3) = a
2664: call VecRestoreArrayF90(x,xx_v,ierr)
2665: .ve
2667: Level: beginner
2669: Note:
2670: If you ONLY intend to read entries from the array and not change any entries you should use `VecGetArrayReadF90()`.
2672: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayReadF90()`
2673: M*/
2675: /*MC
2676: VecGetArrayReadF90 - Accesses a read only array from Fortran. For default PETSc
2677: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2678: this routine is implementation dependent. You MUST call `VecRestoreArrayReadF90()`
2679: when you no longer need access to the array.
2681: Synopsis:
2682: VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2684: Logically Collective
2686: Input Parameter:
2687: . x - vector
2689: Output Parameters:
2690: + xx_v - the Fortran pointer to the array
2691: - ierr - error code
2693: Example of Usage:
2694: .vb
2695: #include <petsc/finclude/petscvec.h>
2696: use petscvec
2698: PetscScalar, pointer :: xx_v(:)
2699: ....
2700: call VecGetArrayReadF90(x,xx_v,ierr)
2701: a = xx_v(3)
2702: call VecRestoreArrayReadF90(x,xx_v,ierr)
2703: .ve
2705: Level: beginner
2707: Note:
2708: If you intend to write entries into the array you must use `VecGetArrayF90()`.
2710: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecGetArrayF90()`
2711: M*/
2713: /*MC
2714: VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2715: `VecGetArrayReadF90()`.
2717: Synopsis:
2718: VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2720: Logically Collective
2722: Input Parameters:
2723: + x - vector
2724: - xx_v - the Fortran pointer to the array
2726: Output Parameter:
2727: . ierr - error code
2729: Example of Usage:
2730: .vb
2731: #include <petsc/finclude/petscvec.h>
2732: use petscvec
2734: PetscScalar, pointer :: xx_v(:)
2735: ....
2736: call VecGetArrayReadF90(x,xx_v,ierr)
2737: a = xx_v(3)
2738: call VecRestoreArrayReadF90(x,xx_v,ierr)
2739: .ve
2741: Level: beginner
2743: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecRestoreArrayF90()`
2744: M*/
2746: /*@C
2747: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2748: processor's portion of the vector data. You MUST call `VecRestoreArray2d()`
2749: when you no longer need access to the array.
2751: Logically Collective
2753: Input Parameters:
2754: + x - the vector
2755: . m - first dimension of two dimensional array
2756: . n - second dimension of two dimensional array
2757: . mstart - first index you will use in first coordinate direction (often 0)
2758: - nstart - first index in the second coordinate direction (often 0)
2760: Output Parameter:
2761: . a - location to put pointer to the array
2763: Level: developer
2765: Notes:
2766: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2767: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2768: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2769: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2771: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2773: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2774: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2775: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2776: @*/
2777: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2778: {
2779: PetscInt i, N;
2780: PetscScalar *aa;
2782: PetscFunctionBegin;
2784: PetscAssertPointer(a, 6);
2786: PetscCall(VecGetLocalSize(x, &N));
2787: 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);
2788: PetscCall(VecGetArray(x, &aa));
2790: PetscCall(PetscMalloc1(m, a));
2791: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2792: *a -= mstart;
2793: PetscFunctionReturn(PETSC_SUCCESS);
2794: }
2796: /*@C
2797: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2798: processor's portion of the vector data. You MUST call `VecRestoreArray2dWrite()`
2799: when you no longer need access to the array.
2801: Logically Collective
2803: Input Parameters:
2804: + x - the vector
2805: . m - first dimension of two dimensional array
2806: . n - second dimension of two dimensional array
2807: . mstart - first index you will use in first coordinate direction (often 0)
2808: - nstart - first index in the second coordinate direction (often 0)
2810: Output Parameter:
2811: . a - location to put pointer to the array
2813: Level: developer
2815: Notes:
2816: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2817: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2818: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2819: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2821: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2823: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2824: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2825: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2826: @*/
2827: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2828: {
2829: PetscInt i, N;
2830: PetscScalar *aa;
2832: PetscFunctionBegin;
2834: PetscAssertPointer(a, 6);
2836: PetscCall(VecGetLocalSize(x, &N));
2837: 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);
2838: PetscCall(VecGetArrayWrite(x, &aa));
2840: PetscCall(PetscMalloc1(m, a));
2841: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2842: *a -= mstart;
2843: PetscFunctionReturn(PETSC_SUCCESS);
2844: }
2846: /*@C
2847: VecRestoreArray2d - Restores a vector after `VecGetArray2d()` has been called.
2849: Logically Collective
2851: Input Parameters:
2852: + x - the vector
2853: . m - first dimension of two dimensional array
2854: . n - second dimension of the two dimensional array
2855: . mstart - first index you will use in first coordinate direction (often 0)
2856: . nstart - first index in the second coordinate direction (often 0)
2857: - a - location of pointer to array obtained from `VecGetArray2d()`
2859: Level: developer
2861: Notes:
2862: For regular PETSc vectors this routine does not involve any copies. For
2863: any special vectors that do not store local vector data in a contiguous
2864: array, this routine will copy the data back into the underlying
2865: vector data structure from the array obtained with `VecGetArray()`.
2867: This routine actually zeros out the `a` pointer.
2869: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2870: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2871: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2872: @*/
2873: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2874: {
2875: void *dummy;
2877: PetscFunctionBegin;
2879: PetscAssertPointer(a, 6);
2881: dummy = (void *)(*a + mstart);
2882: PetscCall(PetscFree(dummy));
2883: PetscCall(VecRestoreArray(x, NULL));
2884: *a = NULL;
2885: PetscFunctionReturn(PETSC_SUCCESS);
2886: }
2888: /*@C
2889: VecRestoreArray2dWrite - Restores a vector after `VecGetArray2dWrite()` has been called.
2891: Logically Collective
2893: Input Parameters:
2894: + x - the vector
2895: . m - first dimension of two dimensional array
2896: . n - second dimension of the two dimensional array
2897: . mstart - first index you will use in first coordinate direction (often 0)
2898: . nstart - first index in the second coordinate direction (often 0)
2899: - a - location of pointer to array obtained from `VecGetArray2d()`
2901: Level: developer
2903: Notes:
2904: For regular PETSc vectors this routine does not involve any copies. For
2905: any special vectors that do not store local vector data in a contiguous
2906: array, this routine will copy the data back into the underlying
2907: vector data structure from the array obtained with `VecGetArray()`.
2909: This routine actually zeros out the `a` pointer.
2911: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2912: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2913: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2914: @*/
2915: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2916: {
2917: void *dummy;
2919: PetscFunctionBegin;
2921: PetscAssertPointer(a, 6);
2923: dummy = (void *)(*a + mstart);
2924: PetscCall(PetscFree(dummy));
2925: PetscCall(VecRestoreArrayWrite(x, NULL));
2926: PetscFunctionReturn(PETSC_SUCCESS);
2927: }
2929: /*@C
2930: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2931: processor's portion of the vector data. You MUST call `VecRestoreArray1d()`
2932: when you no longer need access to the array.
2934: Logically Collective
2936: Input Parameters:
2937: + x - the vector
2938: . m - first dimension of two dimensional array
2939: - mstart - first index you will use in first coordinate direction (often 0)
2941: Output Parameter:
2942: . a - location to put pointer to the array
2944: Level: developer
2946: Notes:
2947: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2948: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2949: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2951: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2953: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2954: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2955: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2956: @*/
2957: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2958: {
2959: PetscInt N;
2961: PetscFunctionBegin;
2963: PetscAssertPointer(a, 4);
2965: PetscCall(VecGetLocalSize(x, &N));
2966: 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);
2967: PetscCall(VecGetArray(x, a));
2968: *a -= mstart;
2969: PetscFunctionReturn(PETSC_SUCCESS);
2970: }
2972: /*@C
2973: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
2974: processor's portion of the vector data. You MUST call `VecRestoreArray1dWrite()`
2975: when you no longer need access to the array.
2977: Logically Collective
2979: Input Parameters:
2980: + x - the vector
2981: . m - first dimension of two dimensional array
2982: - mstart - first index you will use in first coordinate direction (often 0)
2984: Output Parameter:
2985: . a - location to put pointer to the array
2987: Level: developer
2989: Notes:
2990: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2991: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2992: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2994: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2996: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2997: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2998: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2999: @*/
3000: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3001: {
3002: PetscInt N;
3004: PetscFunctionBegin;
3006: PetscAssertPointer(a, 4);
3008: PetscCall(VecGetLocalSize(x, &N));
3009: 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);
3010: PetscCall(VecGetArrayWrite(x, a));
3011: *a -= mstart;
3012: PetscFunctionReturn(PETSC_SUCCESS);
3013: }
3015: /*@C
3016: VecRestoreArray1d - Restores a vector after `VecGetArray1d()` has been called.
3018: Logically Collective
3020: Input Parameters:
3021: + x - the vector
3022: . m - first dimension of two dimensional array
3023: . mstart - first index you will use in first coordinate direction (often 0)
3024: - a - location of pointer to array obtained from `VecGetArray1d()`
3026: Level: developer
3028: Notes:
3029: For regular PETSc vectors this routine does not involve any copies. For
3030: any special vectors that do not store local vector data in a contiguous
3031: array, this routine will copy the data back into the underlying
3032: vector data structure from the array obtained with `VecGetArray1d()`.
3034: This routine actually zeros out the `a` pointer.
3036: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3037: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3038: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3039: @*/
3040: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3041: {
3042: PetscFunctionBegin;
3045: PetscCall(VecRestoreArray(x, NULL));
3046: *a = NULL;
3047: PetscFunctionReturn(PETSC_SUCCESS);
3048: }
3050: /*@C
3051: VecRestoreArray1dWrite - Restores a vector after `VecGetArray1dWrite()` has been called.
3053: Logically Collective
3055: Input Parameters:
3056: + x - the vector
3057: . m - first dimension of two dimensional array
3058: . mstart - first index you will use in first coordinate direction (often 0)
3059: - a - location of pointer to array obtained from `VecGetArray1d()`
3061: Level: developer
3063: Notes:
3064: For regular PETSc vectors this routine does not involve any copies. For
3065: any special vectors that do not store local vector data in a contiguous
3066: array, this routine will copy the data back into the underlying
3067: vector data structure from the array obtained with `VecGetArray1d()`.
3069: This routine actually zeros out the `a` pointer.
3071: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3072: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3073: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3074: @*/
3075: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3076: {
3077: PetscFunctionBegin;
3080: PetscCall(VecRestoreArrayWrite(x, NULL));
3081: *a = NULL;
3082: PetscFunctionReturn(PETSC_SUCCESS);
3083: }
3085: /*@C
3086: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
3087: processor's portion of the vector data. You MUST call `VecRestoreArray3d()`
3088: when you no longer need access to the array.
3090: Logically Collective
3092: Input Parameters:
3093: + x - the vector
3094: . m - first dimension of three dimensional array
3095: . n - second dimension of three dimensional array
3096: . p - third dimension of three dimensional array
3097: . mstart - first index you will use in first coordinate direction (often 0)
3098: . nstart - first index in the second coordinate direction (often 0)
3099: - pstart - first index in the third coordinate direction (often 0)
3101: Output Parameter:
3102: . a - location to put pointer to the array
3104: Level: developer
3106: Notes:
3107: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3108: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3109: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3110: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3112: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3114: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3115: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
3116: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3117: @*/
3118: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3119: {
3120: PetscInt i, N, j;
3121: PetscScalar *aa, **b;
3123: PetscFunctionBegin;
3125: PetscAssertPointer(a, 8);
3127: PetscCall(VecGetLocalSize(x, &N));
3128: 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);
3129: PetscCall(VecGetArray(x, &aa));
3131: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3132: b = (PetscScalar **)((*a) + m);
3133: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3134: for (i = 0; i < m; i++)
3135: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3136: *a -= mstart;
3137: PetscFunctionReturn(PETSC_SUCCESS);
3138: }
3140: /*@C
3141: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
3142: processor's portion of the vector data. You MUST call `VecRestoreArray3dWrite()`
3143: when you no longer need access to the array.
3145: Logically Collective
3147: Input Parameters:
3148: + x - the vector
3149: . m - first dimension of three dimensional array
3150: . n - second dimension of three dimensional array
3151: . p - third dimension of three dimensional array
3152: . mstart - first index you will use in first coordinate direction (often 0)
3153: . nstart - first index in the second coordinate direction (often 0)
3154: - pstart - first index in the third coordinate direction (often 0)
3156: Output Parameter:
3157: . a - location to put pointer to the array
3159: Level: developer
3161: Notes:
3162: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3163: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3164: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3165: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3167: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3169: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3170: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3171: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3172: @*/
3173: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3174: {
3175: PetscInt i, N, j;
3176: PetscScalar *aa, **b;
3178: PetscFunctionBegin;
3180: PetscAssertPointer(a, 8);
3182: PetscCall(VecGetLocalSize(x, &N));
3183: 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);
3184: PetscCall(VecGetArrayWrite(x, &aa));
3186: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3187: b = (PetscScalar **)((*a) + m);
3188: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3189: for (i = 0; i < m; i++)
3190: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3192: *a -= mstart;
3193: PetscFunctionReturn(PETSC_SUCCESS);
3194: }
3196: /*@C
3197: VecRestoreArray3d - Restores a vector after `VecGetArray3d()` has been called.
3199: Logically Collective
3201: Input Parameters:
3202: + x - the vector
3203: . m - first dimension of three dimensional array
3204: . n - second dimension of the three dimensional array
3205: . p - third dimension of the three dimensional array
3206: . mstart - first index you will use in first coordinate direction (often 0)
3207: . nstart - first index in the second coordinate direction (often 0)
3208: . pstart - first index in the third coordinate direction (often 0)
3209: - a - location of pointer to array obtained from VecGetArray3d()
3211: Level: developer
3213: Notes:
3214: For regular PETSc vectors this routine does not involve any copies. For
3215: any special vectors that do not store local vector data in a contiguous
3216: array, this routine will copy the data back into the underlying
3217: vector data structure from the array obtained with `VecGetArray()`.
3219: This routine actually zeros out the `a` pointer.
3221: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3222: `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3223: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3224: @*/
3225: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3226: {
3227: void *dummy;
3229: PetscFunctionBegin;
3231: PetscAssertPointer(a, 8);
3233: dummy = (void *)(*a + mstart);
3234: PetscCall(PetscFree(dummy));
3235: PetscCall(VecRestoreArray(x, NULL));
3236: *a = NULL;
3237: PetscFunctionReturn(PETSC_SUCCESS);
3238: }
3240: /*@C
3241: VecRestoreArray3dWrite - Restores a vector after `VecGetArray3dWrite()` has been called.
3243: Logically Collective
3245: Input Parameters:
3246: + x - the vector
3247: . m - first dimension of three dimensional array
3248: . n - second dimension of the three dimensional array
3249: . p - third dimension of the three dimensional array
3250: . mstart - first index you will use in first coordinate direction (often 0)
3251: . nstart - first index in the second coordinate direction (often 0)
3252: . pstart - first index in the third coordinate direction (often 0)
3253: - a - location of pointer to array obtained from VecGetArray3d()
3255: Level: developer
3257: Notes:
3258: For regular PETSc vectors this routine does not involve any copies. For
3259: any special vectors that do not store local vector data in a contiguous
3260: array, this routine will copy the data back into the underlying
3261: vector data structure from the array obtained with `VecGetArray()`.
3263: This routine actually zeros out the `a` pointer.
3265: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3266: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3267: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3268: @*/
3269: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3270: {
3271: void *dummy;
3273: PetscFunctionBegin;
3275: PetscAssertPointer(a, 8);
3277: dummy = (void *)(*a + mstart);
3278: PetscCall(PetscFree(dummy));
3279: PetscCall(VecRestoreArrayWrite(x, NULL));
3280: *a = NULL;
3281: PetscFunctionReturn(PETSC_SUCCESS);
3282: }
3284: /*@C
3285: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3286: processor's portion of the vector data. You MUST call `VecRestoreArray4d()`
3287: when you no longer need access to the array.
3289: Logically Collective
3291: Input Parameters:
3292: + x - the vector
3293: . m - first dimension of four dimensional array
3294: . n - second dimension of four dimensional array
3295: . p - third dimension of four dimensional array
3296: . q - fourth dimension of four dimensional array
3297: . mstart - first index you will use in first coordinate direction (often 0)
3298: . nstart - first index in the second coordinate direction (often 0)
3299: . pstart - first index in the third coordinate direction (often 0)
3300: - qstart - first index in the fourth coordinate direction (often 0)
3302: Output Parameter:
3303: . a - location to put pointer to the array
3305: Level: developer
3307: Notes:
3308: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3309: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3310: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3311: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3313: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3315: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3316: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3317: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3318: @*/
3319: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3320: {
3321: PetscInt i, N, j, k;
3322: PetscScalar *aa, ***b, **c;
3324: PetscFunctionBegin;
3326: PetscAssertPointer(a, 10);
3328: PetscCall(VecGetLocalSize(x, &N));
3329: 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);
3330: PetscCall(VecGetArray(x, &aa));
3332: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3333: b = (PetscScalar ***)((*a) + m);
3334: c = (PetscScalar **)(b + m * n);
3335: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3336: for (i = 0; i < m; i++)
3337: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3338: for (i = 0; i < m; i++)
3339: for (j = 0; j < n; j++)
3340: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3341: *a -= mstart;
3342: PetscFunctionReturn(PETSC_SUCCESS);
3343: }
3345: /*@C
3346: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3347: processor's portion of the vector data. You MUST call `VecRestoreArray4dWrite()`
3348: when you no longer need access to the array.
3350: Logically Collective
3352: Input Parameters:
3353: + x - the vector
3354: . m - first dimension of four dimensional array
3355: . n - second dimension of four dimensional array
3356: . p - third dimension of four dimensional array
3357: . q - fourth dimension of four dimensional array
3358: . mstart - first index you will use in first coordinate direction (often 0)
3359: . nstart - first index in the second coordinate direction (often 0)
3360: . pstart - first index in the third coordinate direction (often 0)
3361: - qstart - first index in the fourth coordinate direction (often 0)
3363: Output Parameter:
3364: . a - location to put pointer to the array
3366: Level: developer
3368: Notes:
3369: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3370: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3371: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3372: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3374: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3376: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3377: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3378: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3379: @*/
3380: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3381: {
3382: PetscInt i, N, j, k;
3383: PetscScalar *aa, ***b, **c;
3385: PetscFunctionBegin;
3387: PetscAssertPointer(a, 10);
3389: PetscCall(VecGetLocalSize(x, &N));
3390: 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);
3391: PetscCall(VecGetArrayWrite(x, &aa));
3393: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3394: b = (PetscScalar ***)((*a) + m);
3395: c = (PetscScalar **)(b + m * n);
3396: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3397: for (i = 0; i < m; i++)
3398: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3399: for (i = 0; i < m; i++)
3400: for (j = 0; j < n; j++)
3401: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3402: *a -= mstart;
3403: PetscFunctionReturn(PETSC_SUCCESS);
3404: }
3406: /*@C
3407: VecRestoreArray4d - Restores a vector after `VecGetArray4d()` has been called.
3409: Logically Collective
3411: Input Parameters:
3412: + x - the vector
3413: . m - first dimension of four dimensional array
3414: . n - second dimension of the four dimensional array
3415: . p - third dimension of the four dimensional array
3416: . q - fourth dimension of the four dimensional array
3417: . mstart - first index you will use in first coordinate direction (often 0)
3418: . nstart - first index in the second coordinate direction (often 0)
3419: . pstart - first index in the third coordinate direction (often 0)
3420: . qstart - first index in the fourth coordinate direction (often 0)
3421: - a - location of pointer to array obtained from VecGetArray4d()
3423: Level: developer
3425: Notes:
3426: For regular PETSc vectors this routine does not involve any copies. For
3427: any special vectors that do not store local vector data in a contiguous
3428: array, this routine will copy the data back into the underlying
3429: vector data structure from the array obtained with `VecGetArray()`.
3431: This routine actually zeros out the `a` pointer.
3433: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3434: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3435: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecGet`
3436: @*/
3437: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3438: {
3439: void *dummy;
3441: PetscFunctionBegin;
3443: PetscAssertPointer(a, 10);
3445: dummy = (void *)(*a + mstart);
3446: PetscCall(PetscFree(dummy));
3447: PetscCall(VecRestoreArray(x, NULL));
3448: *a = NULL;
3449: PetscFunctionReturn(PETSC_SUCCESS);
3450: }
3452: /*@C
3453: VecRestoreArray4dWrite - Restores a vector after `VecGetArray4dWrite()` has been called.
3455: Logically Collective
3457: Input Parameters:
3458: + x - the vector
3459: . m - first dimension of four dimensional array
3460: . n - second dimension of the four dimensional array
3461: . p - third dimension of the four dimensional array
3462: . q - fourth dimension of the four dimensional array
3463: . mstart - first index you will use in first coordinate direction (often 0)
3464: . nstart - first index in the second coordinate direction (often 0)
3465: . pstart - first index in the third coordinate direction (often 0)
3466: . qstart - first index in the fourth coordinate direction (often 0)
3467: - a - location of pointer to array obtained from `VecGetArray4d()`
3469: Level: developer
3471: Notes:
3472: For regular PETSc vectors this routine does not involve any copies. For
3473: any special vectors that do not store local vector data in a contiguous
3474: array, this routine will copy the data back into the underlying
3475: vector data structure from the array obtained with `VecGetArray()`.
3477: This routine actually zeros out the `a` pointer.
3479: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3480: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3481: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3482: @*/
3483: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3484: {
3485: void *dummy;
3487: PetscFunctionBegin;
3489: PetscAssertPointer(a, 10);
3491: dummy = (void *)(*a + mstart);
3492: PetscCall(PetscFree(dummy));
3493: PetscCall(VecRestoreArrayWrite(x, NULL));
3494: *a = NULL;
3495: PetscFunctionReturn(PETSC_SUCCESS);
3496: }
3498: /*@C
3499: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3500: processor's portion of the vector data. You MUST call `VecRestoreArray2dRead()`
3501: when you no longer need access to the array.
3503: Logically Collective
3505: Input Parameters:
3506: + x - the vector
3507: . m - first dimension of two dimensional array
3508: . n - second dimension of two dimensional array
3509: . mstart - first index you will use in first coordinate direction (often 0)
3510: - nstart - first index in the second coordinate direction (often 0)
3512: Output Parameter:
3513: . a - location to put pointer to the array
3515: Level: developer
3517: Notes:
3518: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
3519: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3520: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3521: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
3523: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3525: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3526: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3527: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3528: @*/
3529: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3530: {
3531: PetscInt i, N;
3532: const PetscScalar *aa;
3534: PetscFunctionBegin;
3536: PetscAssertPointer(a, 6);
3538: PetscCall(VecGetLocalSize(x, &N));
3539: 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);
3540: PetscCall(VecGetArrayRead(x, &aa));
3542: PetscCall(PetscMalloc1(m, a));
3543: for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3544: *a -= mstart;
3545: PetscFunctionReturn(PETSC_SUCCESS);
3546: }
3548: /*@C
3549: VecRestoreArray2dRead - Restores a vector after `VecGetArray2dRead()` has been called.
3551: Logically Collective
3553: Input Parameters:
3554: + x - the vector
3555: . m - first dimension of two dimensional array
3556: . n - second dimension of the two dimensional array
3557: . mstart - first index you will use in first coordinate direction (often 0)
3558: . nstart - first index in the second coordinate direction (often 0)
3559: - a - location of pointer to array obtained from VecGetArray2d()
3561: Level: developer
3563: Notes:
3564: For regular PETSc vectors this routine does not involve any copies. For
3565: any special vectors that do not store local vector data in a contiguous
3566: array, this routine will copy the data back into the underlying
3567: vector data structure from the array obtained with `VecGetArray()`.
3569: This routine actually zeros out the `a` pointer.
3571: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3572: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3573: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3574: @*/
3575: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3576: {
3577: void *dummy;
3579: PetscFunctionBegin;
3581: PetscAssertPointer(a, 6);
3583: dummy = (void *)(*a + mstart);
3584: PetscCall(PetscFree(dummy));
3585: PetscCall(VecRestoreArrayRead(x, NULL));
3586: *a = NULL;
3587: PetscFunctionReturn(PETSC_SUCCESS);
3588: }
3590: /*@C
3591: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3592: processor's portion of the vector data. You MUST call `VecRestoreArray1dRead()`
3593: when you no longer need access to the array.
3595: Logically Collective
3597: Input Parameters:
3598: + x - the vector
3599: . m - first dimension of two dimensional array
3600: - mstart - first index you will use in first coordinate direction (often 0)
3602: Output Parameter:
3603: . a - location to put pointer to the array
3605: Level: developer
3607: Notes:
3608: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3609: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3610: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3612: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3614: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3615: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3616: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3617: @*/
3618: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3619: {
3620: PetscInt N;
3622: PetscFunctionBegin;
3624: PetscAssertPointer(a, 4);
3626: PetscCall(VecGetLocalSize(x, &N));
3627: 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);
3628: PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3629: *a -= mstart;
3630: PetscFunctionReturn(PETSC_SUCCESS);
3631: }
3633: /*@C
3634: VecRestoreArray1dRead - Restores a vector after `VecGetArray1dRead()` has been called.
3636: Logically Collective
3638: Input Parameters:
3639: + x - the vector
3640: . m - first dimension of two dimensional array
3641: . mstart - first index you will use in first coordinate direction (often 0)
3642: - a - location of pointer to array obtained from `VecGetArray1dRead()`
3644: Level: developer
3646: Notes:
3647: For regular PETSc vectors this routine does not involve any copies. For
3648: any special vectors that do not store local vector data in a contiguous
3649: array, this routine will copy the data back into the underlying
3650: vector data structure from the array obtained with `VecGetArray1dRead()`.
3652: This routine actually zeros out the `a` pointer.
3654: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3655: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3656: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3657: @*/
3658: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3659: {
3660: PetscFunctionBegin;
3663: PetscCall(VecRestoreArrayRead(x, NULL));
3664: *a = NULL;
3665: PetscFunctionReturn(PETSC_SUCCESS);
3666: }
3668: /*@C
3669: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3670: processor's portion of the vector data. You MUST call `VecRestoreArray3dRead()`
3671: when you no longer need access to the array.
3673: Logically Collective
3675: Input Parameters:
3676: + x - the vector
3677: . m - first dimension of three dimensional array
3678: . n - second dimension of three dimensional array
3679: . p - third dimension of three dimensional array
3680: . mstart - first index you will use in first coordinate direction (often 0)
3681: . nstart - first index in the second coordinate direction (often 0)
3682: - pstart - first index in the third coordinate direction (often 0)
3684: Output Parameter:
3685: . a - location to put pointer to the array
3687: Level: developer
3689: Notes:
3690: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3691: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3692: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3693: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3dRead()`.
3695: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3697: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3698: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3699: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3700: @*/
3701: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3702: {
3703: PetscInt i, N, j;
3704: const PetscScalar *aa;
3705: PetscScalar **b;
3707: PetscFunctionBegin;
3709: PetscAssertPointer(a, 8);
3711: PetscCall(VecGetLocalSize(x, &N));
3712: 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);
3713: PetscCall(VecGetArrayRead(x, &aa));
3715: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3716: b = (PetscScalar **)((*a) + m);
3717: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3718: for (i = 0; i < m; i++)
3719: for (j = 0; j < n; j++) b[i * n + j] = (PetscScalar *)aa + i * n * p + j * p - pstart;
3720: *a -= mstart;
3721: PetscFunctionReturn(PETSC_SUCCESS);
3722: }
3724: /*@C
3725: VecRestoreArray3dRead - Restores a vector after `VecGetArray3dRead()` has been called.
3727: Logically Collective
3729: Input Parameters:
3730: + x - the vector
3731: . m - first dimension of three dimensional array
3732: . n - second dimension of the three dimensional array
3733: . p - third dimension of the three dimensional array
3734: . mstart - first index you will use in first coordinate direction (often 0)
3735: . nstart - first index in the second coordinate direction (often 0)
3736: . pstart - first index in the third coordinate direction (often 0)
3737: - a - location of pointer to array obtained from `VecGetArray3dRead()`
3739: Level: developer
3741: Notes:
3742: For regular PETSc vectors this routine does not involve any copies. For
3743: any special vectors that do not store local vector data in a contiguous
3744: array, this routine will copy the data back into the underlying
3745: vector data structure from the array obtained with `VecGetArray()`.
3747: This routine actually zeros out the `a` pointer.
3749: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3750: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3751: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3752: @*/
3753: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3754: {
3755: void *dummy;
3757: PetscFunctionBegin;
3759: PetscAssertPointer(a, 8);
3761: dummy = (void *)(*a + mstart);
3762: PetscCall(PetscFree(dummy));
3763: PetscCall(VecRestoreArrayRead(x, NULL));
3764: *a = NULL;
3765: PetscFunctionReturn(PETSC_SUCCESS);
3766: }
3768: /*@C
3769: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3770: processor's portion of the vector data. You MUST call `VecRestoreArray4dRead()`
3771: when you no longer need access to the array.
3773: Logically Collective
3775: Input Parameters:
3776: + x - the vector
3777: . m - first dimension of four dimensional array
3778: . n - second dimension of four dimensional array
3779: . p - third dimension of four dimensional array
3780: . q - fourth dimension of four dimensional array
3781: . mstart - first index you will use in first coordinate direction (often 0)
3782: . nstart - first index in the second coordinate direction (often 0)
3783: . pstart - first index in the third coordinate direction (often 0)
3784: - qstart - first index in the fourth coordinate direction (often 0)
3786: Output Parameter:
3787: . a - location to put pointer to the array
3789: Level: beginner
3791: Notes:
3792: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3793: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3794: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3795: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3797: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3799: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3800: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3801: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3802: @*/
3803: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3804: {
3805: PetscInt i, N, j, k;
3806: const PetscScalar *aa;
3807: PetscScalar ***b, **c;
3809: PetscFunctionBegin;
3811: PetscAssertPointer(a, 10);
3813: PetscCall(VecGetLocalSize(x, &N));
3814: 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);
3815: PetscCall(VecGetArrayRead(x, &aa));
3817: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3818: b = (PetscScalar ***)((*a) + m);
3819: c = (PetscScalar **)(b + m * n);
3820: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3821: for (i = 0; i < m; i++)
3822: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3823: for (i = 0; i < m; i++)
3824: for (j = 0; j < n; j++)
3825: 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;
3826: *a -= mstart;
3827: PetscFunctionReturn(PETSC_SUCCESS);
3828: }
3830: /*@C
3831: VecRestoreArray4dRead - Restores a vector after `VecGetArray4d()` has been called.
3833: Logically Collective
3835: Input Parameters:
3836: + x - the vector
3837: . m - first dimension of four dimensional array
3838: . n - second dimension of the four dimensional array
3839: . p - third dimension of the four dimensional array
3840: . q - fourth dimension of the four dimensional array
3841: . mstart - first index you will use in first coordinate direction (often 0)
3842: . nstart - first index in the second coordinate direction (often 0)
3843: . pstart - first index in the third coordinate direction (often 0)
3844: . qstart - first index in the fourth coordinate direction (often 0)
3845: - a - location of pointer to array obtained from `VecGetArray4dRead()`
3847: Level: beginner
3849: Notes:
3850: For regular PETSc vectors this routine does not involve any copies. For
3851: any special vectors that do not store local vector data in a contiguous
3852: array, this routine will copy the data back into the underlying
3853: vector data structure from the array obtained with `VecGetArray()`.
3855: This routine actually zeros out the `a` pointer.
3857: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3858: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3859: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3860: @*/
3861: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3862: {
3863: void *dummy;
3865: PetscFunctionBegin;
3867: PetscAssertPointer(a, 10);
3869: dummy = (void *)(*a + mstart);
3870: PetscCall(PetscFree(dummy));
3871: PetscCall(VecRestoreArrayRead(x, NULL));
3872: *a = NULL;
3873: PetscFunctionReturn(PETSC_SUCCESS);
3874: }
3876: #if defined(PETSC_USE_DEBUG)
3878: /*@
3879: VecLockGet - Gets the current lock status of a vector
3881: Logically Collective
3883: Input Parameter:
3884: . x - the vector
3886: Output Parameter:
3887: . state - greater than zero indicates the vector is locked for read; less then zero indicates the vector is
3888: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3890: Level: advanced
3892: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3893: @*/
3894: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3895: {
3896: PetscFunctionBegin;
3898: *state = x->lock;
3899: PetscFunctionReturn(PETSC_SUCCESS);
3900: }
3902: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3903: {
3904: PetscFunctionBegin;
3906: PetscAssertPointer(file, 2);
3907: PetscAssertPointer(func, 3);
3908: PetscAssertPointer(line, 4);
3909: #if !PetscDefined(HAVE_THREADSAFETY)
3910: {
3911: const int index = x->lockstack.currentsize - 1;
3913: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupted vec lock stack, have negative index %d", index);
3914: *file = x->lockstack.file[index];
3915: *func = x->lockstack.function[index];
3916: *line = x->lockstack.line[index];
3917: }
3918: #else
3919: *file = NULL;
3920: *func = NULL;
3921: *line = 0;
3922: #endif
3923: PetscFunctionReturn(PETSC_SUCCESS);
3924: }
3926: /*@
3927: VecLockReadPush - Pushes a read-only lock on a vector to prevent it from being written to
3929: Logically Collective
3931: Input Parameter:
3932: . x - the vector
3934: Level: intermediate
3936: Notes:
3937: If this is set then calls to `VecGetArray()` or `VecSetValues()` or any other routines that change the vectors values will generate an error.
3939: The call can be nested, i.e., called multiple times on the same vector, but each `VecLockReadPush()` has to have one matching
3940: `VecLockReadPop()`, which removes the latest read-only lock.
3942: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3943: @*/
3944: PetscErrorCode VecLockReadPush(Vec x)
3945: {
3946: PetscFunctionBegin;
3948: 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");
3949: #if !PetscDefined(HAVE_THREADSAFETY)
3950: {
3951: const char *file, *func;
3952: int index, line;
3954: if ((index = petscstack.currentsize - 2) == -1) {
3955: // vec was locked "outside" of petsc, either in user-land or main. the error message will
3956: // now show this function as the culprit, but it will include the stacktrace
3957: file = "unknown user-file";
3958: func = "unknown_user_function";
3959: line = 0;
3960: } else {
3961: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected petscstack, have negative index %d", index);
3962: file = petscstack.file[index];
3963: func = petscstack.function[index];
3964: line = petscstack.line[index];
3965: }
3966: PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
3967: }
3968: #endif
3969: PetscFunctionReturn(PETSC_SUCCESS);
3970: }
3972: /*@
3973: VecLockReadPop - Pops a read-only lock from a vector
3975: Logically Collective
3977: Input Parameter:
3978: . x - the vector
3980: Level: intermediate
3982: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
3983: @*/
3984: PetscErrorCode VecLockReadPop(Vec x)
3985: {
3986: PetscFunctionBegin;
3988: PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
3989: #if !PetscDefined(HAVE_THREADSAFETY)
3990: {
3991: const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];
3993: PetscStackPop_Private(x->lockstack, previous);
3994: }
3995: #endif
3996: PetscFunctionReturn(PETSC_SUCCESS);
3997: }
3999: /*@C
4000: VecLockWriteSet - Lock or unlock a vector for exclusive read/write access
4002: Logically Collective
4004: Input Parameters:
4005: + x - the vector
4006: - flg - `PETSC_TRUE` to lock the vector for exclusive read/write access; `PETSC_FALSE` to unlock it.
4008: Level: intermediate
4010: Notes:
4011: The function is useful in split-phase computations, which usually have a begin phase and an end phase.
4012: One can call `VecLockWriteSet`(x,`PETSC_TRUE`) in the begin phase to lock a vector for exclusive
4013: access, and call `VecLockWriteSet`(x,`PETSC_FALSE`) in the end phase to unlock the vector from exclusive
4014: access. In this way, one is ensured no other operations can access the vector in between. The code may like
4016: .vb
4017: VecGetArray(x,&xdata); // begin phase
4018: VecLockWriteSet(v,PETSC_TRUE);
4020: Other operations, which can not access x anymore (they can access xdata, of course)
4022: VecRestoreArray(x,&vdata); // end phase
4023: VecLockWriteSet(v,PETSC_FALSE);
4024: .ve
4026: The call can not be nested on the same vector, in other words, one can not call `VecLockWriteSet`(x,`PETSC_TRUE`)
4027: again before calling `VecLockWriteSet`(v,`PETSC_FALSE`).
4029: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
4030: @*/
4031: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
4032: {
4033: PetscFunctionBegin;
4035: if (flg) {
4036: 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");
4037: 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");
4038: x->lock = -1;
4039: } else {
4040: 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");
4041: x->lock = 0;
4042: }
4043: PetscFunctionReturn(PETSC_SUCCESS);
4044: }
4046: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
4047: /*@
4048: VecLockPush - Pushes a read-only lock on a vector to prevent it from being written to
4050: Level: deprecated
4052: .seealso: [](ch_vectors), `Vec`, `VecLockReadPush()`
4053: @*/
4054: PetscErrorCode VecLockPush(Vec x)
4055: {
4056: PetscFunctionBegin;
4057: PetscCall(VecLockReadPush(x));
4058: PetscFunctionReturn(PETSC_SUCCESS);
4059: }
4061: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
4062: /*@
4063: VecLockPop - Pops a read-only lock from a vector
4065: Level: deprecated
4067: .seealso: [](ch_vectors), `Vec`, `VecLockReadPop()`
4068: @*/
4069: PetscErrorCode VecLockPop(Vec x)
4070: {
4071: PetscFunctionBegin;
4072: PetscCall(VecLockReadPop(x));
4073: PetscFunctionReturn(PETSC_SUCCESS);
4074: }
4076: #endif