Actual source code: rvector.c
1: /*
2: Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
3: These are the vector functions the user calls.
4: */
5: #include "petsc/private/sfimpl.h"
6: #include "petscsystypes.h"
7: #include <petsc/private/vecimpl.h>
9: PetscInt VecGetSubVectorSavedStateId = -1;
11: #if PetscDefined(USE_DEBUG)
12: // this is a no-op '0' macro in optimized builds
13: PetscErrorCode VecValidValues_Internal(Vec vec, PetscInt argnum, PetscBool begin)
14: {
15: PetscFunctionBegin;
16: if (vec->petscnative || vec->ops->getarray) {
17: PetscInt n;
18: const PetscScalar *x;
19: PetscOffloadMask mask;
21: PetscCall(VecGetOffloadMask(vec, &mask));
22: if (!PetscOffloadHost(mask)) PetscFunctionReturn(PETSC_SUCCESS);
23: PetscCall(VecGetLocalSize(vec, &n));
24: PetscCall(VecGetArrayRead(vec, &x));
25: for (PetscInt i = 0; i < n; i++) {
26: if (begin) {
27: PetscCheck(!PetscIsInfOrNanScalar(x[i]), PETSC_COMM_SELF, PETSC_ERR_FP, "Vec entry at local location %" PetscInt_FMT " is not-a-number or infinite at beginning of function: Parameter number %" PetscInt_FMT, i, argnum);
28: } else {
29: PetscCheck(!PetscIsInfOrNanScalar(x[i]), PETSC_COMM_SELF, PETSC_ERR_FP, "Vec entry at local location %" PetscInt_FMT " is not-a-number or infinite at end of function: Parameter number %" PetscInt_FMT, i, argnum);
30: }
31: }
32: PetscCall(VecRestoreArrayRead(vec, &x));
33: }
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
36: #endif
38: /*@
39: VecMaxPointwiseDivide - Computes the maximum of the componentwise division `max = max_i abs(x[i]/y[i])`.
41: Logically Collective
43: Input Parameters:
44: + x - the numerators
45: - y - the denominators
47: Output Parameter:
48: . max - the result
50: Level: advanced
52: Notes:
53: `x` and `y` may be the same vector
55: if a particular `y[i]` is zero, it is treated as 1 in the above formula
57: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`
58: @*/
59: PetscErrorCode VecMaxPointwiseDivide(Vec x, Vec y, PetscReal *max)
60: {
61: PetscFunctionBegin;
64: PetscAssertPointer(max, 3);
67: PetscCheckSameTypeAndComm(x, 1, y, 2);
68: VecCheckSameSize(x, 1, y, 2);
69: VecCheckAssembled(x);
70: VecCheckAssembled(y);
71: PetscCall(VecLockReadPush(x));
72: PetscCall(VecLockReadPush(y));
73: PetscUseTypeMethod(x, maxpointwisedivide, y, max);
74: PetscCall(VecLockReadPop(x));
75: PetscCall(VecLockReadPop(y));
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: /*@
80: VecDot - Computes the vector dot product.
82: Collective
84: Input Parameters:
85: + x - first vector
86: - y - second vector
88: Output Parameter:
89: . val - the dot product
91: Level: intermediate
93: Notes for Users of Complex Numbers:
94: For complex vectors, `VecDot()` computes
95: $ val = (x,y) = y^H x,
96: where y^H denotes the conjugate transpose of y. Note that this corresponds to the usual "mathematicians" complex
97: inner product where the SECOND argument gets the complex conjugate. Since the `BLASdot()` complex conjugates the first
98: first argument we call the `BLASdot()` with the arguments reversed.
100: Use `VecTDot()` for the indefinite form
101: $ val = (x,y) = y^T x,
102: where y^T denotes the transpose of y.
104: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDotRealPart()`
105: @*/
106: PetscErrorCode VecDot(Vec x, Vec y, PetscScalar *val)
107: {
108: PetscFunctionBegin;
111: PetscAssertPointer(val, 3);
114: PetscCheckSameTypeAndComm(x, 1, y, 2);
115: VecCheckSameSize(x, 1, y, 2);
116: VecCheckAssembled(x);
117: VecCheckAssembled(y);
119: PetscCall(VecLockReadPush(x));
120: PetscCall(VecLockReadPush(y));
121: PetscCall(PetscLogEventBegin(VEC_Dot, x, y, 0, 0));
122: PetscUseTypeMethod(x, dot, y, val);
123: PetscCall(PetscLogEventEnd(VEC_Dot, x, y, 0, 0));
124: PetscCall(VecLockReadPop(x));
125: PetscCall(VecLockReadPop(y));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: /*@
130: VecDotRealPart - Computes the real part of the vector dot product.
132: Collective
134: Input Parameters:
135: + x - first vector
136: - y - second vector
138: Output Parameter:
139: . val - the real part of the dot product;
141: Level: intermediate
143: Notes for Users of Complex Numbers:
144: See `VecDot()` for more details on the definition of the dot product for complex numbers
146: For real numbers this returns the same value as `VecDot()`
148: For complex numbers in C^n (that is a vector of n components with a complex number for each component) this is equal to the usual real dot product on the
149: the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)
151: Developer Notes:
152: This is not currently optimized to compute only the real part of the dot product.
154: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDot()`, `VecDotNorm2()`
155: @*/
156: PetscErrorCode VecDotRealPart(Vec x, Vec y, PetscReal *val)
157: {
158: PetscScalar fdot;
160: PetscFunctionBegin;
161: PetscCall(VecDot(x, y, &fdot));
162: *val = PetscRealPart(fdot);
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: /*@
167: VecNorm - Computes the vector norm.
169: Collective
171: Input Parameters:
172: + x - the vector
173: - type - the type of the norm requested
175: Output Parameter:
176: . val - the norm
178: Level: intermediate
180: Notes:
181: See `NormType` for descriptions of each norm.
183: For complex numbers `NORM_1` will return the traditional 1 norm of the 2 norm of the complex
184: numbers; that is the 1 norm of the absolute values of the complex entries. In PETSc 3.6 and
185: earlier releases it returned the 1 norm of the 1 norm of the complex entries (what is
186: returned by the BLAS routine `asum()`). Both are valid norms but most people expect the former.
188: This routine stashes the computed norm value, repeated calls before the vector entries are
189: changed are then rapid since the precomputed value is immediately available. Certain vector
190: operations such as `VecSet()` store the norms so the value is immediately available and does
191: not need to be explicitly computed. `VecScale()` updates any stashed norm values, thus calls
192: after `VecScale()` do not need to explicitly recompute the norm.
194: .seealso: [](ch_vectors), `Vec`, `NormType`, `VecDot()`, `VecTDot()`, `VecDotBegin()`, `VecDotEnd()`, `VecNormAvailable()`,
195: `VecNormBegin()`, `VecNormEnd()`, `NormType()`
196: @*/
197: PetscErrorCode VecNorm(Vec x, NormType type, PetscReal *val)
198: {
199: PetscBool flg = PETSC_TRUE;
201: PetscFunctionBegin;
204: VecCheckAssembled(x);
206: PetscAssertPointer(val, 3);
208: PetscCall(VecNormAvailable(x, type, &flg, val));
209: // check that all MPI processes call this routine together and have same availability
210: if (PetscDefined(USE_DEBUG)) {
211: PetscMPIInt b0 = (PetscMPIInt)flg, b1[2], b2[2];
212: b1[0] = -b0;
213: b1[1] = b0;
214: PetscCall(MPIU_Allreduce(b1, b2, 2, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
215: PetscCheck(-b2[0] == b2[1], PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_WRONGSTATE, "Some MPI processes have cached norm, others do not. This may happen when some MPI processes call VecGetArray() and some others do not.");
216: }
217: if (flg) PetscFunctionReturn(PETSC_SUCCESS);
219: PetscCall(VecLockReadPush(x));
220: PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
221: PetscUseTypeMethod(x, norm, type, val);
222: PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
223: PetscCall(VecLockReadPop(x));
225: if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: /*@
230: VecNormAvailable - Returns the vector norm if it is already known. That is, it has been previously computed and cached in the vector
232: Not Collective
234: Input Parameters:
235: + x - the vector
236: - 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
237: `NORM_1_AND_2`, which computes both norms and stores them
238: in a two element array.
240: Output Parameters:
241: + available - `PETSC_TRUE` if the val returned is valid
242: - val - the norm
244: Level: intermediate
246: Developer Notes:
247: `PETSC_HAVE_SLOW_BLAS_NORM2` will cause a C (loop unrolled) version of the norm to be used, rather
248: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
249: (as opposed to vendor provided) because the FORTRAN BLAS `NRM2()` routine is very slow.
251: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
252: `VecNormBegin()`, `VecNormEnd()`
253: @*/
254: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
255: {
256: PetscFunctionBegin;
259: PetscAssertPointer(available, 3);
260: PetscAssertPointer(val, 4);
262: if (type == NORM_1_AND_2) {
263: *available = PETSC_FALSE;
264: } else {
265: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
266: }
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: /*@
271: VecNormalize - Normalizes a vector by its 2-norm.
273: Collective
275: Input Parameter:
276: . x - the vector
278: Output Parameter:
279: . val - the vector norm before normalization. May be `NULL` if the value is not needed.
281: Level: intermediate
283: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
284: @*/
285: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
286: {
287: PetscReal norm;
289: PetscFunctionBegin;
292: PetscCall(VecSetErrorIfLocked(x, 1));
293: if (val) PetscAssertPointer(val, 2);
294: PetscCall(PetscLogEventBegin(VEC_Normalize, x, 0, 0, 0));
295: PetscCall(VecNorm(x, NORM_2, &norm));
296: if (norm == 0.0) PetscCall(PetscInfo(x, "Vector of zero norm can not be normalized; Returning only the zero norm\n"));
297: else if (PetscIsInfOrNanReal(norm)) PetscCall(PetscInfo(x, "Vector with Inf or Nan norm can not be normalized; Returning only the norm\n"));
298: else {
299: PetscScalar s = 1.0 / norm;
300: PetscCall(VecScale(x, s));
301: }
302: PetscCall(PetscLogEventEnd(VEC_Normalize, x, 0, 0, 0));
303: if (val) *val = norm;
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }
307: /*@C
308: VecMax - Determines the vector component with maximum real part and its location.
310: Collective
312: Input Parameter:
313: . x - the vector
315: Output Parameters:
316: + p - the index of `val` (pass `NULL` if you don't want this) in the vector
317: - val - the maximum component
319: Level: intermediate
321: Notes:
322: Returns the value `PETSC_MIN_REAL` and negative `p` if the vector is of length 0.
324: Returns the smallest index with the maximum value
326: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
327: @*/
328: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
329: {
330: PetscFunctionBegin;
333: VecCheckAssembled(x);
334: if (p) PetscAssertPointer(p, 2);
335: PetscAssertPointer(val, 3);
336: PetscCall(VecLockReadPush(x));
337: PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
338: PetscUseTypeMethod(x, max, p, val);
339: PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
340: PetscCall(VecLockReadPop(x));
341: PetscFunctionReturn(PETSC_SUCCESS);
342: }
344: /*@C
345: VecMin - Determines the vector component with minimum real part and its location.
347: Collective
349: Input Parameter:
350: . x - the vector
352: Output Parameters:
353: + p - the index of `val` (pass `NULL` if you don't want this location) in the vector
354: - val - the minimum component
356: Level: intermediate
358: Notes:
359: Returns the value `PETSC_MAX_REAL` and negative `p` if the vector is of length 0.
361: This returns the smallest index with the minimum value
363: .seealso: [](ch_vectors), `Vec`, `VecMax()`
364: @*/
365: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
366: {
367: PetscFunctionBegin;
370: VecCheckAssembled(x);
371: if (p) PetscAssertPointer(p, 2);
372: PetscAssertPointer(val, 3);
373: PetscCall(VecLockReadPush(x));
374: PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
375: PetscUseTypeMethod(x, min, p, val);
376: PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
377: PetscCall(VecLockReadPop(x));
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: /*@
382: VecTDot - Computes an indefinite vector dot product. That is, this
383: routine does NOT use the complex conjugate.
385: Collective
387: Input Parameters:
388: + x - first vector
389: - y - second vector
391: Output Parameter:
392: . val - the dot product
394: Level: intermediate
396: Notes for Users of Complex Numbers:
397: For complex vectors, `VecTDot()` computes the indefinite form
398: $ val = (x,y) = y^T x,
399: where y^T denotes the transpose of y.
401: Use `VecDot()` for the inner product
402: $ val = (x,y) = y^H x,
403: where y^H denotes the conjugate transpose of y.
405: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
406: @*/
407: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
408: {
409: PetscFunctionBegin;
412: PetscAssertPointer(val, 3);
415: PetscCheckSameTypeAndComm(x, 1, y, 2);
416: VecCheckSameSize(x, 1, y, 2);
417: VecCheckAssembled(x);
418: VecCheckAssembled(y);
420: PetscCall(VecLockReadPush(x));
421: PetscCall(VecLockReadPush(y));
422: PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
423: PetscUseTypeMethod(x, tdot, y, val);
424: PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
425: PetscCall(VecLockReadPop(x));
426: PetscCall(VecLockReadPop(y));
427: PetscFunctionReturn(PETSC_SUCCESS);
428: }
430: PetscErrorCode VecScaleAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
431: {
432: PetscReal norms[4];
433: PetscBool flgs[4];
434: PetscScalar one = 1.0;
436: PetscFunctionBegin;
439: VecCheckAssembled(x);
440: PetscCall(VecSetErrorIfLocked(x, 1));
441: if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);
443: /* get current stashed norms */
444: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
446: PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
447: VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
448: PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));
450: PetscCall(PetscObjectStateIncrease((PetscObject)x));
451: /* put the scaled stashed norms back into the Vec */
452: for (PetscInt i = 0; i < 4; i++) {
453: PetscReal ar = PetscAbsScalar(alpha);
454: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
455: }
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }
459: /*@
460: VecScale - Scales a vector.
462: Not Collective
464: Input Parameters:
465: + x - the vector
466: - alpha - the scalar
468: Level: intermediate
470: Note:
471: For a vector with n components, `VecScale()` computes x[i] = alpha * x[i], for i=1,...,n.
473: .seealso: [](ch_vectors), `Vec`, `VecSet()`
474: @*/
475: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
476: {
477: PetscFunctionBegin;
478: PetscCall(VecScaleAsync_Private(x, alpha, NULL));
479: PetscFunctionReturn(PETSC_SUCCESS);
480: }
482: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
483: {
484: PetscFunctionBegin;
487: VecCheckAssembled(x);
489: PetscCall(VecSetErrorIfLocked(x, 1));
491: if (alpha == 0) {
492: PetscReal norm;
493: PetscBool set;
495: PetscCall(VecNormAvailable(x, NORM_2, &set, &norm));
496: if (set == PETSC_TRUE && norm == 0) PetscFunctionReturn(PETSC_SUCCESS);
497: }
498: PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
499: VecMethodDispatch(x, dctx, VecAsyncFnName(Set), set, (Vec, PetscScalar, PetscDeviceContext), alpha);
500: PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
501: PetscCall(PetscObjectStateIncrease((PetscObject)x));
503: /* norms can be simply set (if |alpha|*N not too large) */
505: {
506: PetscReal val = PetscAbsScalar(alpha);
507: const PetscInt N = x->map->N;
509: if (N == 0) {
510: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0));
511: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
512: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
513: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
514: } else if (val > PETSC_MAX_REAL / N) {
515: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
516: } else {
517: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
518: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
519: val *= PetscSqrtReal((PetscReal)N);
520: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
521: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
522: }
523: }
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: /*@
528: VecSet - Sets all components of a vector to a single scalar value.
530: Logically Collective
532: Input Parameters:
533: + x - the vector
534: - alpha - the scalar
536: Level: beginner
538: Notes:
539: For a vector of dimension n, `VecSet()` sets x[i] = alpha, for i=1,...,n,
540: so that all vector entries then equal the identical
541: scalar value, `alpha`. Use the more general routine
542: `VecSetValues()` to set different vector entries.
544: You CANNOT call this after you have called `VecSetValues()` but before you call
545: `VecAssemblyBegin()`
547: If `alpha` is zero and the norm of the vector is known to be zero then this skips the unneeded zeroing process
549: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
550: @*/
551: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
552: {
553: PetscFunctionBegin;
554: PetscCall(VecSetAsync_Private(x, alpha, NULL));
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }
558: PetscErrorCode VecAXPYAsync_Private(Vec y, PetscScalar alpha, Vec x, PetscDeviceContext dctx)
559: {
560: PetscFunctionBegin;
565: PetscCheckSameTypeAndComm(x, 3, y, 1);
566: VecCheckSameSize(x, 3, y, 1);
567: VecCheckAssembled(x);
568: VecCheckAssembled(y);
570: if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
571: PetscCall(VecSetErrorIfLocked(y, 1));
572: if (x == y) {
573: PetscCall(VecScale(y, alpha + 1.0));
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
576: PetscCall(VecLockReadPush(x));
577: PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
578: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPY), axpy, (Vec, PetscScalar, Vec, PetscDeviceContext), alpha, x);
579: PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
580: PetscCall(VecLockReadPop(x));
581: PetscCall(PetscObjectStateIncrease((PetscObject)y));
582: PetscFunctionReturn(PETSC_SUCCESS);
583: }
584: /*@
585: VecAXPY - Computes `y = alpha x + y`.
587: Logically Collective
589: Input Parameters:
590: + alpha - the scalar
591: . x - vector scale by `alpha`
592: - y - vector accumulated into
594: Output Parameter:
595: . y - output vector
597: Level: intermediate
599: Notes:
600: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
601: .vb
602: VecAXPY(y,alpha,x) y = alpha x + y
603: VecAYPX(y,beta,x) y = x + beta y
604: VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
605: VecWAXPY(w,alpha,x,y) w = alpha x + y
606: VecAXPBYPCZ(z,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
607: VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
608: .ve
610: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
611: @*/
612: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
613: {
614: PetscFunctionBegin;
615: PetscCall(VecAXPYAsync_Private(y, alpha, x, NULL));
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: PetscErrorCode VecAYPXAsync_Private(Vec y, PetscScalar beta, Vec x, PetscDeviceContext dctx)
620: {
621: PetscFunctionBegin;
626: PetscCheckSameTypeAndComm(x, 3, y, 1);
627: VecCheckSameSize(x, 1, y, 3);
628: VecCheckAssembled(x);
629: VecCheckAssembled(y);
631: PetscCall(VecSetErrorIfLocked(y, 1));
632: if (x == y) {
633: PetscCall(VecScale(y, beta + 1.0));
634: PetscFunctionReturn(PETSC_SUCCESS);
635: }
636: PetscCall(VecLockReadPush(x));
637: if (beta == (PetscScalar)0.0) {
638: PetscCall(VecCopy(x, y));
639: } else {
640: PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
641: VecMethodDispatch(y, dctx, VecAsyncFnName(AYPX), aypx, (Vec, PetscScalar, Vec, PetscDeviceContext), beta, x);
642: PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
643: PetscCall(PetscObjectStateIncrease((PetscObject)y));
644: }
645: PetscCall(VecLockReadPop(x));
646: PetscFunctionReturn(PETSC_SUCCESS);
647: }
649: /*@
650: VecAYPX - Computes `y = x + beta y`.
652: Logically Collective
654: Input Parameters:
655: + beta - the scalar
656: . x - the unscaled vector
657: - y - the vector to be scaled
659: Output Parameter:
660: . y - output vector
662: Level: intermediate
664: Developer Notes:
665: The implementation is optimized for `beta` of -1.0, 0.0, and 1.0
667: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
668: @*/
669: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
670: {
671: PetscFunctionBegin;
672: PetscCall(VecAYPXAsync_Private(y, beta, x, NULL));
673: PetscFunctionReturn(PETSC_SUCCESS);
674: }
676: PetscErrorCode VecAXPBYAsync_Private(Vec y, PetscScalar alpha, PetscScalar beta, Vec x, PetscDeviceContext dctx)
677: {
678: PetscFunctionBegin;
683: PetscCheckSameTypeAndComm(x, 4, y, 1);
684: VecCheckSameSize(y, 1, x, 4);
685: VecCheckAssembled(x);
686: VecCheckAssembled(y);
689: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
690: if (x == y) {
691: PetscCall(VecScale(y, alpha + beta));
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: PetscCall(VecSetErrorIfLocked(y, 1));
696: PetscCall(VecLockReadPush(x));
697: PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
698: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPBY), axpby, (Vec, PetscScalar, PetscScalar, Vec, PetscDeviceContext), alpha, beta, x);
699: PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
700: PetscCall(PetscObjectStateIncrease((PetscObject)y));
701: PetscCall(VecLockReadPop(x));
702: PetscFunctionReturn(PETSC_SUCCESS);
703: }
704: /*@
705: VecAXPBY - Computes `y = alpha x + beta y`.
707: Logically Collective
709: Input Parameters:
710: + alpha - first scalar
711: . beta - second scalar
712: . x - the first scaled vector
713: - y - the second scaled vector
715: Output Parameter:
716: . y - output vector
718: Level: intermediate
720: Developer Notes:
721: The implementation is optimized for `alpha` and/or `beta` values of 0.0 and 1.0
723: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
724: @*/
725: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
726: {
727: PetscFunctionBegin;
728: PetscCall(VecAXPBYAsync_Private(y, alpha, beta, x, NULL));
729: PetscFunctionReturn(PETSC_SUCCESS);
730: }
732: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
733: {
734: PetscFunctionBegin;
741: PetscCheckSameTypeAndComm(x, 5, y, 6);
742: PetscCheckSameTypeAndComm(x, 5, z, 1);
743: VecCheckSameSize(x, 5, y, 6);
744: VecCheckSameSize(x, 5, z, 1);
745: PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
746: PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
747: VecCheckAssembled(x);
748: VecCheckAssembled(y);
749: VecCheckAssembled(z);
753: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
755: PetscCall(VecSetErrorIfLocked(z, 1));
756: PetscCall(VecLockReadPush(x));
757: PetscCall(VecLockReadPush(y));
758: PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
759: VecMethodDispatch(z, dctx, VecAsyncFnName(AXPBYPCZ), axpbypcz, (Vec, PetscScalar, PetscScalar, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, beta, gamma, x, y);
760: PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
761: PetscCall(PetscObjectStateIncrease((PetscObject)z));
762: PetscCall(VecLockReadPop(x));
763: PetscCall(VecLockReadPop(y));
764: PetscFunctionReturn(PETSC_SUCCESS);
765: }
766: /*@
767: VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`
769: Logically Collective
771: Input Parameters:
772: + alpha - first scalar
773: . beta - second scalar
774: . gamma - third scalar
775: . x - first vector
776: . y - second vector
777: - z - third vector
779: Output Parameter:
780: . z - output vector
782: Level: intermediate
784: Note:
785: `x`, `y` and `z` must be different vectors
787: Developer Notes:
788: The implementation is optimized for `alpha` of 1.0 and `gamma` of 1.0 or 0.0
790: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
791: @*/
792: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
793: {
794: PetscFunctionBegin;
795: PetscCall(VecAXPBYPCZAsync_Private(z, alpha, beta, gamma, x, y, NULL));
796: PetscFunctionReturn(PETSC_SUCCESS);
797: }
799: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
800: {
801: PetscFunctionBegin;
808: PetscCheckSameTypeAndComm(x, 3, y, 4);
809: PetscCheckSameTypeAndComm(y, 4, w, 1);
810: VecCheckSameSize(x, 3, y, 4);
811: VecCheckSameSize(x, 3, w, 1);
812: PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
813: PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
814: VecCheckAssembled(x);
815: VecCheckAssembled(y);
817: PetscCall(VecSetErrorIfLocked(w, 1));
819: PetscCall(VecLockReadPush(x));
820: PetscCall(VecLockReadPush(y));
821: if (alpha == (PetscScalar)0.0) {
822: PetscCall(VecCopyAsync_Private(y, w, dctx));
823: } else {
824: PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
825: VecMethodDispatch(w, dctx, VecAsyncFnName(WAXPY), waxpy, (Vec, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, x, y);
826: PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
827: PetscCall(PetscObjectStateIncrease((PetscObject)w));
828: }
829: PetscCall(VecLockReadPop(x));
830: PetscCall(VecLockReadPop(y));
831: PetscFunctionReturn(PETSC_SUCCESS);
832: }
834: /*@
835: VecWAXPY - Computes `w = alpha x + y`.
837: Logically Collective
839: Input Parameters:
840: + alpha - the scalar
841: . x - first vector, multiplied by `alpha`
842: - y - second vector
844: Output Parameter:
845: . w - the result
847: Level: intermediate
849: Note:
850: `w` cannot be either `x` or `y`, but `x` and `y` can be the same
852: Developer Notes:
853: The implementation is optimized for alpha of -1.0, 0.0, and 1.0
855: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
856: @*/
857: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
858: {
859: PetscFunctionBegin;
860: PetscCall(VecWAXPYAsync_Private(w, alpha, x, y, NULL));
861: PetscFunctionReturn(PETSC_SUCCESS);
862: }
864: /*@C
865: VecSetValues - Inserts or adds values into certain locations of a vector.
867: Not Collective
869: Input Parameters:
870: + x - vector to insert in
871: . ni - number of elements to add
872: . ix - indices where to add
873: . y - array of values
874: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries
876: Level: beginner
878: Notes:
879: .vb
880: `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
881: .ve
883: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
884: options cannot be mixed without intervening calls to the assembly
885: routines.
887: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
888: MUST be called after all calls to `VecSetValues()` have been completed.
890: VecSetValues() uses 0-based indices in Fortran as well as in C.
892: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
893: negative indices may be passed in ix. These rows are
894: simply ignored. This allows easily inserting element load matrices
895: with homogeneous Dirichlet boundary conditions that you don't want represented
896: in the vector.
898: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
899: `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
900: @*/
901: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
902: {
903: PetscFunctionBeginHot;
905: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
906: PetscAssertPointer(ix, 3);
907: PetscAssertPointer(y, 4);
910: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
911: PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
912: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
913: PetscCall(PetscObjectStateIncrease((PetscObject)x));
914: PetscFunctionReturn(PETSC_SUCCESS);
915: }
917: /*@C
918: VecGetValues - Gets values from certain locations of a vector. Currently
919: can only get values on the same processor on which they are owned
921: Not Collective
923: Input Parameters:
924: + x - vector to get values from
925: . ni - number of elements to get
926: - ix - indices where to get them from (in global 1d numbering)
928: Output Parameter:
929: . y - array of values
931: Level: beginner
933: Notes:
934: The user provides the allocated array y; it is NOT allocated in this routine
936: `VecGetValues()` gets y[i] = x[ix[i]], for i=0,...,ni-1.
938: `VecAssemblyBegin()` and `VecAssemblyEnd()` MUST be called before calling this if `VecSetValues()` or related routine has been called
940: VecGetValues() uses 0-based indices in Fortran as well as in C.
942: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
943: negative indices may be passed in ix. These rows are
944: simply ignored.
946: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
947: @*/
948: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
949: {
950: PetscFunctionBegin;
952: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
953: PetscAssertPointer(ix, 3);
954: PetscAssertPointer(y, 4);
956: VecCheckAssembled(x);
957: PetscUseTypeMethod(x, getvalues, ni, ix, y);
958: PetscFunctionReturn(PETSC_SUCCESS);
959: }
961: /*@C
962: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
964: Not Collective
966: Input Parameters:
967: + x - vector to insert in
968: . ni - number of blocks to add
969: . ix - indices where to add in block count, rather than element count
970: . y - array of values
971: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries
973: Level: intermediate
975: Notes:
976: `VecSetValuesBlocked()` sets x[bs*ix[i]+j] = y[bs*i+j],
977: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
979: Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
980: options cannot be mixed without intervening calls to the assembly
981: routines.
983: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
984: MUST be called after all calls to `VecSetValuesBlocked()` have been completed.
986: `VecSetValuesBlocked()` uses 0-based indices in Fortran as well as in C.
988: Negative indices may be passed in ix, these rows are
989: simply ignored. This allows easily inserting element load matrices
990: with homogeneous Dirichlet boundary conditions that you don't want represented
991: in the vector.
993: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
994: `VecSetValues()`
995: @*/
996: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
997: {
998: PetscFunctionBeginHot;
1000: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1001: PetscAssertPointer(ix, 3);
1002: PetscAssertPointer(y, 4);
1005: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1006: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1007: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1008: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1009: PetscFunctionReturn(PETSC_SUCCESS);
1010: }
1012: /*@C
1013: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1014: using a local ordering of the nodes.
1016: Not Collective
1018: Input Parameters:
1019: + x - vector to insert in
1020: . ni - number of elements to add
1021: . ix - indices where to add
1022: . y - array of values
1023: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1025: Level: intermediate
1027: Notes:
1028: `VecSetValuesLocal()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
1030: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
1031: options cannot be mixed without intervening calls to the assembly
1032: routines.
1034: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1035: MUST be called after all calls to `VecSetValuesLocal()` have been completed.
1037: `VecSetValuesLocal()` uses 0-based indices in Fortran as well as in C.
1039: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1040: `VecSetValuesBlockedLocal()`
1041: @*/
1042: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1043: {
1044: PetscInt lixp[128], *lix = lixp;
1046: PetscFunctionBeginHot;
1048: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1049: PetscAssertPointer(ix, 3);
1050: PetscAssertPointer(y, 4);
1053: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1054: if (!x->ops->setvalueslocal) {
1055: if (x->map->mapping) {
1056: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1057: PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1058: PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1059: if (ni > 128) PetscCall(PetscFree(lix));
1060: } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1061: } else PetscUseTypeMethod(x, setvalueslocal, ni, ix, y, iora);
1062: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1063: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1064: PetscFunctionReturn(PETSC_SUCCESS);
1065: }
1067: /*@
1068: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1069: using a local ordering of the nodes.
1071: Not Collective
1073: Input Parameters:
1074: + x - vector to insert in
1075: . ni - number of blocks to add
1076: . ix - indices where to add in block count, not element count
1077: . y - array of values
1078: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1080: Level: intermediate
1082: Notes:
1083: `VecSetValuesBlockedLocal()` sets x[bs*ix[i]+j] = y[bs*i+j],
1084: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with `VecSetBlockSize()`.
1086: Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1087: options cannot be mixed without intervening calls to the assembly
1088: routines.
1090: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1091: MUST be called after all calls to `VecSetValuesBlockedLocal()` have been completed.
1093: `VecSetValuesBlockedLocal()` uses 0-based indices in Fortran as well as in C.
1095: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1096: `VecSetLocalToGlobalMapping()`
1097: @*/
1098: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1099: {
1100: PetscInt lixp[128], *lix = lixp;
1102: PetscFunctionBeginHot;
1104: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1105: PetscAssertPointer(ix, 3);
1106: PetscAssertPointer(y, 4);
1108: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1109: if (x->map->mapping) {
1110: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1111: PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1112: PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1113: if (ni > 128) PetscCall(PetscFree(lix));
1114: } else {
1115: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1116: }
1117: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1118: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1119: PetscFunctionReturn(PETSC_SUCCESS);
1120: }
1122: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1123: {
1124: PetscFunctionBegin;
1127: VecCheckAssembled(x);
1129: if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1130: PetscAssertPointer(y, 3);
1131: for (PetscInt i = 0; i < nv; ++i) {
1134: PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1135: VecCheckSameSize(x, 1, y[i], 3);
1136: VecCheckAssembled(y[i]);
1137: PetscCall(VecLockReadPush(y[i]));
1138: }
1139: PetscAssertPointer(result, 4);
1142: PetscCall(VecLockReadPush(x));
1143: PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1144: PetscCall((*mxdot)(x, nv, y, result));
1145: PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1146: PetscCall(VecLockReadPop(x));
1147: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1148: PetscFunctionReturn(PETSC_SUCCESS);
1149: }
1151: /*@
1152: VecMTDot - Computes indefinite vector multiple dot products.
1153: That is, it does NOT use the complex conjugate.
1155: Collective
1157: Input Parameters:
1158: + x - one vector
1159: . nv - number of vectors
1160: - y - array of vectors. Note that vectors are pointers
1162: Output Parameter:
1163: . val - array of the dot products
1165: Level: intermediate
1167: Notes for Users of Complex Numbers:
1168: For complex vectors, `VecMTDot()` computes the indefinite form
1169: $ val = (x,y) = y^T x,
1170: where y^T denotes the transpose of y.
1172: Use `VecMDot()` for the inner product
1173: $ val = (x,y) = y^H x,
1174: where y^H denotes the conjugate transpose of y.
1176: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1177: @*/
1178: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1179: {
1180: PetscFunctionBegin;
1182: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1183: PetscFunctionReturn(PETSC_SUCCESS);
1184: }
1186: /*@
1187: VecMDot - Computes multiple vector dot products.
1189: Collective
1191: Input Parameters:
1192: + x - one vector
1193: . nv - number of vectors
1194: - y - array of vectors.
1196: Output Parameter:
1197: . val - array of the dot products (does not allocate the array)
1199: Level: intermediate
1201: Notes for Users of Complex Numbers:
1202: For complex vectors, `VecMDot()` computes
1203: $ val = (x,y) = y^H x,
1204: where y^H denotes the conjugate transpose of y.
1206: Use `VecMTDot()` for the indefinite form
1207: $ val = (x,y) = y^T x,
1208: where y^T denotes the transpose of y.
1210: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`
1211: @*/
1212: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1213: {
1214: PetscFunctionBegin;
1216: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1217: PetscFunctionReturn(PETSC_SUCCESS);
1218: }
1220: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
1221: {
1222: PetscFunctionBegin;
1224: VecCheckAssembled(y);
1226: PetscCall(VecSetErrorIfLocked(y, 1));
1227: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1228: if (nv) {
1229: PetscInt zeros = 0;
1231: PetscAssertPointer(alpha, 3);
1232: PetscAssertPointer(x, 4);
1233: for (PetscInt i = 0; i < nv; ++i) {
1237: PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1238: VecCheckSameSize(y, 1, x[i], 4);
1239: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1240: VecCheckAssembled(x[i]);
1241: PetscCall(VecLockReadPush(x[i]));
1242: zeros += alpha[i] == (PetscScalar)0.0;
1243: }
1245: if (zeros < nv) {
1246: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1247: VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1248: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1249: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1250: }
1252: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1253: }
1254: PetscFunctionReturn(PETSC_SUCCESS);
1255: }
1257: /*@
1258: VecMAXPY - Computes `y = y + sum alpha[i] x[i]`
1260: Logically Collective
1262: Input Parameters:
1263: + nv - number of scalars and x-vectors
1264: . alpha - array of scalars
1265: . y - one vector
1266: - x - array of vectors
1268: Level: intermediate
1270: Note:
1271: `y` cannot be any of the `x` vectors
1273: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`,`VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1274: @*/
1275: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1276: {
1277: PetscFunctionBegin;
1278: PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1279: PetscFunctionReturn(PETSC_SUCCESS);
1280: }
1282: /*@
1283: VecMAXPBY - Computes `y = beta y + sum alpha[i] x[i]`
1285: Logically Collective
1287: Input Parameters:
1288: + nv - number of scalars and x-vectors
1289: . alpha - array of scalars
1290: . beta - scalar
1291: . y - one vector
1292: - x - array of vectors
1294: Level: intermediate
1296: Note:
1297: `y` cannot be any of the `x` vectors.
1299: Developer Notes:
1300: This is a convenience routine, but implementations might be able to optimize it, for example, when `beta` is zero.
1302: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1303: @*/
1304: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
1305: {
1306: PetscFunctionBegin;
1308: VecCheckAssembled(y);
1310: PetscCall(VecSetErrorIfLocked(y, 1));
1311: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1314: if (y->ops->maxpby) {
1315: PetscInt zeros = 0;
1317: if (nv) {
1318: PetscAssertPointer(alpha, 3);
1319: PetscAssertPointer(x, 5);
1320: }
1322: for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
1326: PetscCheckSameTypeAndComm(y, 1, x[i], 5);
1327: VecCheckSameSize(y, 1, x[i], 5);
1328: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1329: VecCheckAssembled(x[i]);
1330: PetscCall(VecLockReadPush(x[i]));
1331: zeros += alpha[i] == (PetscScalar)0.0;
1332: }
1334: if (zeros < nv) { // has nonzero alpha
1335: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1336: PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1337: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1338: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1339: } else {
1340: PetscCall(VecScale(y, beta));
1341: }
1343: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1344: } else { // no maxpby
1345: if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1346: else PetscCall(VecScale(y, beta));
1347: PetscCall(VecMAXPY(y, nv, alpha, x));
1348: }
1349: PetscFunctionReturn(PETSC_SUCCESS);
1350: }
1352: /*@
1353: VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
1354: in the order they appear in the array. The concatenated vector resides on the same
1355: communicator and is the same type as the source vectors.
1357: Collective
1359: Input Parameters:
1360: + nx - number of vectors to be concatenated
1361: - X - array containing the vectors to be concatenated in the order of concatenation
1363: Output Parameters:
1364: + Y - concatenated vector
1365: - x_is - array of index sets corresponding to the concatenated components of `Y` (pass `NULL` if not needed)
1367: Level: advanced
1369: Notes:
1370: Concatenation is similar to the functionality of a `VECNEST` object; they both represent combination of
1371: different vector spaces. However, concatenated vectors do not store any information about their
1372: sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1373: manipulation of data in the concatenated vector that corresponds to the original components at creation.
1375: This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1376: has to operate on combined vector spaces and cannot utilize `VECNEST` objects due to incompatibility with
1377: bound projections.
1379: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1380: @*/
1381: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1382: {
1383: MPI_Comm comm;
1384: VecType vec_type;
1385: Vec Ytmp, Xtmp;
1386: IS *is_tmp;
1387: PetscInt i, shift = 0, Xnl, Xng, Xbegin;
1389: PetscFunctionBegin;
1393: PetscAssertPointer(Y, 3);
1395: if ((*X)->ops->concatenate) {
1396: /* use the dedicated concatenation function if available */
1397: PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1398: } else {
1399: /* loop over vectors and start creating IS */
1400: comm = PetscObjectComm((PetscObject)(*X));
1401: PetscCall(VecGetType(*X, &vec_type));
1402: PetscCall(PetscMalloc1(nx, &is_tmp));
1403: for (i = 0; i < nx; i++) {
1404: PetscCall(VecGetSize(X[i], &Xng));
1405: PetscCall(VecGetLocalSize(X[i], &Xnl));
1406: PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1407: PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1408: shift += Xng;
1409: }
1410: /* create the concatenated vector */
1411: PetscCall(VecCreate(comm, &Ytmp));
1412: PetscCall(VecSetType(Ytmp, vec_type));
1413: PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1414: PetscCall(VecSetUp(Ytmp));
1415: /* copy data from X array to Y and return */
1416: for (i = 0; i < nx; i++) {
1417: PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1418: PetscCall(VecCopy(X[i], Xtmp));
1419: PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1420: }
1421: *Y = Ytmp;
1422: if (x_is) {
1423: *x_is = is_tmp;
1424: } else {
1425: for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1426: PetscCall(PetscFree(is_tmp));
1427: }
1428: }
1429: PetscFunctionReturn(PETSC_SUCCESS);
1430: }
1432: /* A helper function for VecGetSubVector to check if we can implement it with no-copy (i.e. the subvector shares
1433: memory with the original vector), and the block size of the subvector.
1435: Input Parameters:
1436: + X - the original vector
1437: - is - the index set of the subvector
1439: Output Parameters:
1440: + contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
1441: . start - start of contiguous block, as an offset from the start of the ownership range of the original vector
1442: - blocksize - the block size of the subvector
1444: */
1445: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1446: {
1447: PetscInt gstart, gend, lstart;
1448: PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1449: PetscInt n, N, ibs, vbs, bs = -1;
1451: PetscFunctionBegin;
1452: PetscCall(ISGetLocalSize(is, &n));
1453: PetscCall(ISGetSize(is, &N));
1454: PetscCall(ISGetBlockSize(is, &ibs));
1455: PetscCall(VecGetBlockSize(X, &vbs));
1456: PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1457: PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1458: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1459: if (ibs > 1) {
1460: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1461: bs = ibs;
1462: } else {
1463: if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1464: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1465: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1466: }
1468: *contig = red[0];
1469: *start = lstart;
1470: *blocksize = bs;
1471: PetscFunctionReturn(PETSC_SUCCESS);
1472: }
1474: /* A helper function for VecGetSubVector, to be used when we have to build a standalone subvector through VecScatter
1476: Input Parameters:
1477: + X - the original vector
1478: . is - the index set of the subvector
1479: - bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()
1481: Output Parameter:
1482: . Z - the subvector, which will compose the VecScatter context on output
1483: */
1484: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1485: {
1486: PetscInt n, N;
1487: VecScatter vscat;
1488: Vec Y;
1490: PetscFunctionBegin;
1491: PetscCall(ISGetLocalSize(is, &n));
1492: PetscCall(ISGetSize(is, &N));
1493: PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1494: PetscCall(VecSetSizes(Y, n, N));
1495: PetscCall(VecSetBlockSize(Y, bs));
1496: PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1497: PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1498: PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1499: PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1500: PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1501: PetscCall(VecScatterDestroy(&vscat));
1502: *Z = Y;
1503: PetscFunctionReturn(PETSC_SUCCESS);
1504: }
1506: /*@
1507: VecGetSubVector - Gets a vector representing part of another vector
1509: Collective
1511: Input Parameters:
1512: + X - vector from which to extract a subvector
1513: - is - index set representing portion of `X` to extract
1515: Output Parameter:
1516: . Y - subvector corresponding to `is`
1518: Level: advanced
1520: Notes:
1521: The subvector `Y` should be returned with `VecRestoreSubVector()`.
1522: `X` and must be defined on the same communicator
1524: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1525: modifying the subvector. Other non-overlapping subvectors can still be obtained from X using this function.
1527: 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`.
1529: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1530: @*/
1531: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1532: {
1533: Vec Z;
1535: PetscFunctionBegin;
1538: PetscCheckSameComm(X, 1, is, 2);
1539: PetscAssertPointer(Y, 3);
1540: if (X->ops->getsubvector) {
1541: PetscUseTypeMethod(X, getsubvector, is, &Z);
1542: } else { /* Default implementation currently does no caching */
1543: PetscBool contig;
1544: PetscInt n, N, start, bs;
1546: PetscCall(ISGetLocalSize(is, &n));
1547: PetscCall(ISGetSize(is, &N));
1548: PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1549: if (contig) { /* We can do a no-copy implementation */
1550: const PetscScalar *x;
1551: PetscInt state = 0;
1552: PetscBool isstd, iscuda, iship;
1554: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1555: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1556: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1557: if (iscuda) {
1558: #if defined(PETSC_HAVE_CUDA)
1559: const PetscScalar *x_d;
1560: PetscMPIInt size;
1561: PetscOffloadMask flg;
1563: PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1564: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1565: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1566: if (x) x += start;
1567: if (x_d) x_d += start;
1568: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1569: if (size == 1) {
1570: PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1571: } else {
1572: PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1573: }
1574: Z->offloadmask = flg;
1575: #endif
1576: } else if (iship) {
1577: #if defined(PETSC_HAVE_HIP)
1578: const PetscScalar *x_d;
1579: PetscMPIInt size;
1580: PetscOffloadMask flg;
1582: PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1583: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1584: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1585: if (x) x += start;
1586: if (x_d) x_d += start;
1587: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1588: if (size == 1) {
1589: PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1590: } else {
1591: PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1592: }
1593: Z->offloadmask = flg;
1594: #endif
1595: } else if (isstd) {
1596: PetscMPIInt size;
1598: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1599: PetscCall(VecGetArrayRead(X, &x));
1600: if (x) x += start;
1601: if (size == 1) {
1602: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1603: } else {
1604: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1605: }
1606: PetscCall(VecRestoreArrayRead(X, &x));
1607: } else { /* default implementation: use place array */
1608: PetscCall(VecGetArrayRead(X, &x));
1609: PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1610: PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1611: PetscCall(VecSetSizes(Z, n, N));
1612: PetscCall(VecSetBlockSize(Z, bs));
1613: PetscCall(VecPlaceArray(Z, x ? x + start : NULL));
1614: PetscCall(VecRestoreArrayRead(X, &x));
1615: }
1617: /* this is relevant only in debug mode */
1618: PetscCall(VecLockGet(X, &state));
1619: if (state) PetscCall(VecLockReadPush(Z));
1620: Z->ops->placearray = NULL;
1621: Z->ops->replacearray = NULL;
1622: } else { /* Have to create a scatter and do a copy */
1623: PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1624: }
1625: }
1626: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1627: if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1628: PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1629: *Y = Z;
1630: PetscFunctionReturn(PETSC_SUCCESS);
1631: }
1633: /*@
1634: VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`
1636: Collective
1638: Input Parameters:
1639: + X - vector from which subvector was obtained
1640: . is - index set representing the subset of `X`
1641: - Y - subvector being restored
1643: Level: advanced
1645: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1646: @*/
1647: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1648: {
1649: PETSC_UNUSED PetscObjectState dummystate = 0;
1650: PetscBool unchanged;
1652: PetscFunctionBegin;
1655: PetscCheckSameComm(X, 1, is, 2);
1656: PetscAssertPointer(Y, 3);
1659: if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1660: else {
1661: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1662: if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1663: VecScatter scatter;
1664: PetscInt state;
1666: PetscCall(VecLockGet(X, &state));
1667: PetscCheck(state == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec X is locked for read-only or read/write access");
1669: PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1670: if (scatter) {
1671: PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1672: PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1673: } else {
1674: PetscBool iscuda, iship;
1675: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1676: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1678: if (iscuda) {
1679: #if defined(PETSC_HAVE_CUDA)
1680: PetscOffloadMask ymask = (*Y)->offloadmask;
1682: /* The offloadmask of X dictates where to move memory
1683: If X GPU data is valid, then move Y data on GPU if needed
1684: Otherwise, move back to the CPU */
1685: switch (X->offloadmask) {
1686: case PETSC_OFFLOAD_BOTH:
1687: if (ymask == PETSC_OFFLOAD_CPU) {
1688: PetscCall(VecCUDAResetArray(*Y));
1689: } else if (ymask == PETSC_OFFLOAD_GPU) {
1690: X->offloadmask = PETSC_OFFLOAD_GPU;
1691: }
1692: break;
1693: case PETSC_OFFLOAD_GPU:
1694: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1695: break;
1696: case PETSC_OFFLOAD_CPU:
1697: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1698: break;
1699: case PETSC_OFFLOAD_UNALLOCATED:
1700: case PETSC_OFFLOAD_KOKKOS:
1701: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1702: }
1703: #endif
1704: } else if (iship) {
1705: #if defined(PETSC_HAVE_HIP)
1706: PetscOffloadMask ymask = (*Y)->offloadmask;
1708: /* The offloadmask of X dictates where to move memory
1709: If X GPU data is valid, then move Y data on GPU if needed
1710: Otherwise, move back to the CPU */
1711: switch (X->offloadmask) {
1712: case PETSC_OFFLOAD_BOTH:
1713: if (ymask == PETSC_OFFLOAD_CPU) {
1714: PetscCall(VecHIPResetArray(*Y));
1715: } else if (ymask == PETSC_OFFLOAD_GPU) {
1716: X->offloadmask = PETSC_OFFLOAD_GPU;
1717: }
1718: break;
1719: case PETSC_OFFLOAD_GPU:
1720: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1721: break;
1722: case PETSC_OFFLOAD_CPU:
1723: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1724: break;
1725: case PETSC_OFFLOAD_UNALLOCATED:
1726: case PETSC_OFFLOAD_KOKKOS:
1727: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1728: }
1729: #endif
1730: } else {
1731: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1732: PetscCall(VecResetArray(*Y));
1733: }
1734: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1735: }
1736: }
1737: }
1738: PetscCall(VecDestroy(Y));
1739: PetscFunctionReturn(PETSC_SUCCESS);
1740: }
1742: /*@
1743: VecCreateLocalVector - Creates a vector object suitable for use with `VecGetLocalVector()` and friends. You must call `VecDestroy()` when the
1744: vector is no longer needed.
1746: Not Collective.
1748: Input Parameter:
1749: . v - The vector for which the local vector is desired.
1751: Output Parameter:
1752: . w - Upon exit this contains the local vector.
1754: Level: beginner
1756: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1757: @*/
1758: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1759: {
1760: PetscMPIInt size;
1762: PetscFunctionBegin;
1764: PetscAssertPointer(w, 2);
1765: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1766: if (size == 1) PetscCall(VecDuplicate(v, w));
1767: else if (v->ops->createlocalvector) PetscUseTypeMethod(v, createlocalvector, w);
1768: else {
1769: VecType type;
1770: PetscInt n;
1772: PetscCall(VecCreate(PETSC_COMM_SELF, w));
1773: PetscCall(VecGetLocalSize(v, &n));
1774: PetscCall(VecSetSizes(*w, n, n));
1775: PetscCall(VecGetBlockSize(v, &n));
1776: PetscCall(VecSetBlockSize(*w, n));
1777: PetscCall(VecGetType(v, &type));
1778: PetscCall(VecSetType(*w, type));
1779: }
1780: PetscFunctionReturn(PETSC_SUCCESS);
1781: }
1783: /*@
1784: VecGetLocalVectorRead - Maps the local portion of a vector into a
1785: vector.
1787: Not Collective.
1789: Input Parameter:
1790: . v - The vector for which the local vector is desired.
1792: Output Parameter:
1793: . w - Upon exit this contains the local vector.
1795: Level: beginner
1797: Notes:
1798: You must call `VecRestoreLocalVectorRead()` when the local
1799: vector is no longer needed.
1801: This function is similar to `VecGetArrayRead()` which maps the local
1802: portion into a raw pointer. `VecGetLocalVectorRead()` is usually
1803: almost as efficient as `VecGetArrayRead()` but in certain circumstances
1804: `VecGetLocalVectorRead()` can be much more efficient than
1805: `VecGetArrayRead()`. This is because the construction of a contiguous
1806: array representing the vector data required by `VecGetArrayRead()` can
1807: be an expensive operation for certain vector types. For example, for
1808: GPU vectors `VecGetArrayRead()` requires that the data between device
1809: and host is synchronized.
1811: Unlike `VecGetLocalVector()`, this routine is not collective and
1812: preserves cached information.
1814: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1815: @*/
1816: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1817: {
1818: PetscFunctionBegin;
1821: VecCheckSameLocalSize(v, 1, w, 2);
1822: if (v->ops->getlocalvectorread) {
1823: PetscUseTypeMethod(v, getlocalvectorread, w);
1824: } else {
1825: PetscScalar *a;
1827: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1828: PetscCall(VecPlaceArray(w, a));
1829: }
1830: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1831: PetscCall(VecLockReadPush(v));
1832: PetscCall(VecLockReadPush(w));
1833: PetscFunctionReturn(PETSC_SUCCESS);
1834: }
1836: /*@
1837: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1838: previously mapped into a vector using `VecGetLocalVectorRead()`.
1840: Not Collective.
1842: Input Parameters:
1843: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
1844: - w - The vector into which the local portion of `v` was mapped.
1846: Level: beginner
1848: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1849: @*/
1850: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1851: {
1852: PetscFunctionBegin;
1855: if (v->ops->restorelocalvectorread) {
1856: PetscUseTypeMethod(v, restorelocalvectorread, w);
1857: } else {
1858: const PetscScalar *a;
1860: PetscCall(VecGetArrayRead(w, &a));
1861: PetscCall(VecRestoreArrayRead(v, &a));
1862: PetscCall(VecResetArray(w));
1863: }
1864: PetscCall(VecLockReadPop(v));
1865: PetscCall(VecLockReadPop(w));
1866: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1867: PetscFunctionReturn(PETSC_SUCCESS);
1868: }
1870: /*@
1871: VecGetLocalVector - Maps the local portion of a vector into a
1872: vector.
1874: Collective
1876: Input Parameter:
1877: . v - The vector for which the local vector is desired.
1879: Output Parameter:
1880: . w - Upon exit this contains the local vector.
1882: Level: beginner
1884: Notes:
1885: You must call `VecRestoreLocalVector()` when the local
1886: vector is no longer needed.
1888: This function is similar to `VecGetArray()` which maps the local
1889: portion into a raw pointer. `VecGetLocalVector()` is usually about as
1890: efficient as `VecGetArray()` but in certain circumstances
1891: `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1892: This is because the construction of a contiguous array representing
1893: the vector data required by `VecGetArray()` can be an expensive
1894: operation for certain vector types. For example, for GPU vectors
1895: `VecGetArray()` requires that the data between device and host is
1896: synchronized.
1898: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1899: @*/
1900: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1901: {
1902: PetscFunctionBegin;
1905: VecCheckSameLocalSize(v, 1, w, 2);
1906: if (v->ops->getlocalvector) {
1907: PetscUseTypeMethod(v, getlocalvector, w);
1908: } else {
1909: PetscScalar *a;
1911: PetscCall(VecGetArray(v, &a));
1912: PetscCall(VecPlaceArray(w, a));
1913: }
1914: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1915: PetscFunctionReturn(PETSC_SUCCESS);
1916: }
1918: /*@
1919: VecRestoreLocalVector - Unmaps the local portion of a vector
1920: previously mapped into a vector using `VecGetLocalVector()`.
1922: Logically Collective.
1924: Input Parameters:
1925: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
1926: - w - The vector into which the local portion of `v` was mapped.
1928: Level: beginner
1930: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1931: @*/
1932: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1933: {
1934: PetscFunctionBegin;
1937: if (v->ops->restorelocalvector) {
1938: PetscUseTypeMethod(v, restorelocalvector, w);
1939: } else {
1940: PetscScalar *a;
1941: PetscCall(VecGetArray(w, &a));
1942: PetscCall(VecRestoreArray(v, &a));
1943: PetscCall(VecResetArray(w));
1944: }
1945: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1946: PetscCall(PetscObjectStateIncrease((PetscObject)v));
1947: PetscFunctionReturn(PETSC_SUCCESS);
1948: }
1950: /*@C
1951: VecGetArray - Returns a pointer to a contiguous array that contains this
1952: MPI processes's portion of the vector data
1954: Logically Collective
1956: Input Parameter:
1957: . x - the vector
1959: Output Parameter:
1960: . a - location to put pointer to the array
1962: Level: beginner
1964: Notes:
1965: For the standard PETSc vectors, `VecGetArray()` returns a pointer to the local data array and
1966: does not use any copies. If the underlying vector data is not stored in a contiguous array
1967: this routine will copy the data to a contiguous array and return a pointer to that. You MUST
1968: call `VecRestoreArray()` when you no longer need access to the array.
1970: Fortran Notes:
1971: `VecGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayF90()`
1973: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
1974: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
1975: @*/
1976: PetscErrorCode VecGetArray(Vec x, PetscScalar **a)
1977: {
1978: PetscFunctionBegin;
1980: PetscCall(VecSetErrorIfLocked(x, 1));
1981: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
1982: PetscUseTypeMethod(x, getarray, a);
1983: } else if (x->petscnative) { /* VECSTANDARD */
1984: *a = *((PetscScalar **)x->data);
1985: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
1986: PetscFunctionReturn(PETSC_SUCCESS);
1987: }
1989: /*@C
1990: VecRestoreArray - Restores a vector after `VecGetArray()` has been called and the array is no longer needed
1992: Logically Collective
1994: Input Parameters:
1995: + x - the vector
1996: - a - location of pointer to array obtained from `VecGetArray()`
1998: Level: beginner
2000: Fortran Notes:
2001: `VecRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayF90()`
2003: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2004: `VecGetArrayPair()`, `VecRestoreArrayPair()`
2005: @*/
2006: PetscErrorCode VecRestoreArray(Vec x, PetscScalar **a)
2007: {
2008: PetscFunctionBegin;
2010: if (a) PetscAssertPointer(a, 2);
2011: if (x->ops->restorearray) {
2012: PetscUseTypeMethod(x, restorearray, a);
2013: } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
2014: if (a) *a = NULL;
2015: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2016: PetscFunctionReturn(PETSC_SUCCESS);
2017: }
2018: /*@C
2019: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
2021: Not Collective
2023: Input Parameter:
2024: . x - the vector
2026: Output Parameter:
2027: . a - the array
2029: Level: beginner
2031: Notes:
2032: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2034: Unlike `VecGetArray()`, preserves cached information like vector norms.
2036: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
2037: implementations may require a copy, but such implementations should cache the contiguous representation so that
2038: only one copy is performed when this routine is called multiple times in sequence.
2040: Fortran Notes:
2041: `VecGetArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayReadF90()`
2043: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2044: @*/
2045: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar **a)
2046: {
2047: PetscFunctionBegin;
2049: PetscAssertPointer(a, 2);
2050: if (x->ops->getarrayread) {
2051: PetscUseTypeMethod(x, getarrayread, a);
2052: } else if (x->ops->getarray) {
2053: PetscObjectState state;
2055: /* VECNEST, VECCUDA, VECKOKKOS etc */
2056: // x->ops->getarray may bump the object state, but since we know this is a read-only get
2057: // we can just undo that
2058: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2059: PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
2060: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2061: } else if (x->petscnative) {
2062: /* VECSTANDARD */
2063: *a = *((PetscScalar **)x->data);
2064: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
2065: PetscFunctionReturn(PETSC_SUCCESS);
2066: }
2068: /*@C
2069: VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`
2071: Not Collective
2073: Input Parameters:
2074: + x - the vector
2075: - a - the array
2077: Level: beginner
2079: Fortran Notes:
2080: `VecRestoreArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayReadF90()`
2082: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2083: @*/
2084: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar **a)
2085: {
2086: PetscFunctionBegin;
2088: if (a) PetscAssertPointer(a, 2);
2089: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
2090: /* nothing */
2091: } else if (x->ops->restorearrayread) { /* VECNEST */
2092: PetscUseTypeMethod(x, restorearrayread, a);
2093: } else { /* No one? */
2094: PetscObjectState state;
2096: // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2097: // we can just undo that
2098: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2099: PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2100: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2101: }
2102: if (a) *a = NULL;
2103: PetscFunctionReturn(PETSC_SUCCESS);
2104: }
2106: /*@C
2107: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
2108: MPI processes's portion of the vector data.
2110: Logically Collective
2112: Input Parameter:
2113: . x - the vector
2115: Output Parameter:
2116: . a - location to put pointer to the array
2118: Level: intermediate
2120: Note:
2121: The values in this array are NOT valid, the caller of this routine is responsible for putting
2122: values into the array; any values it does not set will be invalid.
2124: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2126: For vectors associated with GPUs, the host and device vectors are not synchronized before
2127: giving access. If you need correct values in the array use `VecGetArray()`
2129: Fortran Notes:
2130: `VecGetArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayWriteF90()`
2132: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteF90()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
2133: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
2134: @*/
2135: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar **a)
2136: {
2137: PetscFunctionBegin;
2139: PetscAssertPointer(a, 2);
2140: PetscCall(VecSetErrorIfLocked(x, 1));
2141: if (x->ops->getarraywrite) {
2142: PetscUseTypeMethod(x, getarraywrite, a);
2143: } else {
2144: PetscCall(VecGetArray(x, a));
2145: }
2146: PetscFunctionReturn(PETSC_SUCCESS);
2147: }
2149: /*@C
2150: VecRestoreArrayWrite - Restores a vector after `VecGetArrayWrite()` has been called.
2152: Logically Collective
2154: Input Parameters:
2155: + x - the vector
2156: - a - location of pointer to array obtained from `VecGetArray()`
2158: Level: beginner
2160: Fortran Notes:
2161: `VecRestoreArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayWriteF90()`
2163: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteF90()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2164: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2165: @*/
2166: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar **a)
2167: {
2168: PetscFunctionBegin;
2170: if (a) PetscAssertPointer(a, 2);
2171: if (x->ops->restorearraywrite) {
2172: PetscUseTypeMethod(x, restorearraywrite, a);
2173: } else if (x->ops->restorearray) {
2174: PetscUseTypeMethod(x, restorearray, a);
2175: }
2176: if (a) *a = NULL;
2177: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2178: PetscFunctionReturn(PETSC_SUCCESS);
2179: }
2181: /*@C
2182: VecGetArrays - Returns a pointer to the arrays in a set of vectors
2183: that were created by a call to `VecDuplicateVecs()`.
2185: Logically Collective; No Fortran Support
2187: Input Parameters:
2188: + x - the vectors
2189: - n - the number of vectors
2191: Output Parameter:
2192: . a - location to put pointer to the array
2194: Level: intermediate
2196: Note:
2197: You MUST call `VecRestoreArrays()` when you no longer need access to the arrays.
2199: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2200: @*/
2201: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2202: {
2203: PetscInt i;
2204: PetscScalar **q;
2206: PetscFunctionBegin;
2207: PetscAssertPointer(x, 1);
2209: PetscAssertPointer(a, 3);
2210: PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2211: PetscCall(PetscMalloc1(n, &q));
2212: for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2213: *a = q;
2214: PetscFunctionReturn(PETSC_SUCCESS);
2215: }
2217: /*@C
2218: VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2219: has been called.
2221: Logically Collective; No Fortran Support
2223: Input Parameters:
2224: + x - the vector
2225: . n - the number of vectors
2226: - a - location of pointer to arrays obtained from `VecGetArrays()`
2228: Notes:
2229: For regular PETSc vectors this routine does not involve any copies. For
2230: any special vectors that do not store local vector data in a contiguous
2231: array, this routine will copy the data back into the underlying
2232: vector data structure from the arrays obtained with `VecGetArrays()`.
2234: Level: intermediate
2236: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2237: @*/
2238: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2239: {
2240: PetscInt i;
2241: PetscScalar **q = *a;
2243: PetscFunctionBegin;
2244: PetscAssertPointer(x, 1);
2246: PetscAssertPointer(a, 3);
2248: for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2249: PetscCall(PetscFree(q));
2250: PetscFunctionReturn(PETSC_SUCCESS);
2251: }
2253: /*@C
2254: VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g.,
2255: `VECCUDA`), the returned pointer will be a device pointer to the device memory that contains
2256: this MPI processes's portion of the vector data.
2258: Logically Collective; No Fortran Support
2260: Input Parameter:
2261: . x - the vector
2263: Output Parameters:
2264: + a - location to put pointer to the array
2265: - mtype - memory type of the array
2267: Level: beginner
2269: Note:
2270: Device data is guaranteed to have the latest value. Otherwise, when this is a host vector
2271: (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host
2272: pointer.
2274: For `VECKOKKOS`, if Kokkos is configured without device (e.g., use serial or openmp), per
2275: this function, the vector works like `VECSEQ`/`VECMPI`; otherwise, it works like `VECCUDA` or
2276: `VECHIP` etc.
2278: Use `VecRestoreArrayAndMemType()` when the array access is no longer needed.
2280: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`,
2281: `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2282: @*/
2283: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2284: {
2285: PetscFunctionBegin;
2288: PetscAssertPointer(a, 2);
2289: if (mtype) PetscAssertPointer(mtype, 3);
2290: PetscCall(VecSetErrorIfLocked(x, 1));
2291: if (x->ops->getarrayandmemtype) {
2292: /* VECCUDA, VECKOKKOS etc */
2293: PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2294: } else {
2295: /* VECSTANDARD, VECNEST, VECVIENNACL */
2296: PetscCall(VecGetArray(x, a));
2297: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2298: }
2299: PetscFunctionReturn(PETSC_SUCCESS);
2300: }
2302: /*@C
2303: VecRestoreArrayAndMemType - Restores a vector after `VecGetArrayAndMemType()` has been called.
2305: Logically Collective; No Fortran Support
2307: Input Parameters:
2308: + x - the vector
2309: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`
2311: Level: beginner
2313: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`,
2314: `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2315: @*/
2316: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar **a)
2317: {
2318: PetscFunctionBegin;
2321: if (a) PetscAssertPointer(a, 2);
2322: if (x->ops->restorearrayandmemtype) {
2323: /* VECCUDA, VECKOKKOS etc */
2324: PetscUseTypeMethod(x, restorearrayandmemtype, a);
2325: } else {
2326: /* VECNEST, VECVIENNACL */
2327: PetscCall(VecRestoreArray(x, a));
2328: } /* VECSTANDARD does nothing */
2329: if (a) *a = NULL;
2330: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2331: PetscFunctionReturn(PETSC_SUCCESS);
2332: }
2334: /*@C
2335: VecGetArrayReadAndMemType - Like `VecGetArrayRead()`, but if the input vector is a device vector, it will return a read-only device pointer.
2336: The returned pointer is guaranteed to point to up-to-date data. For host vectors, it functions as `VecGetArrayRead()`.
2338: Not Collective; No Fortran Support
2340: Input Parameter:
2341: . x - the vector
2343: Output Parameters:
2344: + a - the array
2345: - mtype - memory type of the array
2347: Level: beginner
2349: Notes:
2350: The array must be returned using a matching call to `VecRestoreArrayReadAndMemType()`.
2352: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2353: @*/
2354: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar **a, PetscMemType *mtype)
2355: {
2356: PetscFunctionBegin;
2359: PetscAssertPointer(a, 2);
2360: if (mtype) PetscAssertPointer(mtype, 3);
2361: if (x->ops->getarrayreadandmemtype) {
2362: /* VECCUDA/VECHIP though they are also petscnative */
2363: PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2364: } else if (x->ops->getarrayandmemtype) {
2365: /* VECKOKKOS */
2366: PetscObjectState state;
2368: // see VecGetArrayRead() for why
2369: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2370: PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2371: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2372: } else {
2373: PetscCall(VecGetArrayRead(x, a));
2374: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2375: }
2376: PetscFunctionReturn(PETSC_SUCCESS);
2377: }
2379: /*@C
2380: VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`
2382: Not Collective; No Fortran Support
2384: Input Parameters:
2385: + x - the vector
2386: - a - the array
2388: Level: beginner
2390: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2391: @*/
2392: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar **a)
2393: {
2394: PetscFunctionBegin;
2397: if (a) PetscAssertPointer(a, 2);
2398: if (x->ops->restorearrayreadandmemtype) {
2399: /* VECCUDA/VECHIP */
2400: PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2401: } else if (!x->petscnative) {
2402: /* VECNEST */
2403: PetscCall(VecRestoreArrayRead(x, a));
2404: }
2405: if (a) *a = NULL;
2406: PetscFunctionReturn(PETSC_SUCCESS);
2407: }
2409: /*@C
2410: VecGetArrayWriteAndMemType - Like `VecGetArrayWrite()`, but if this is a device vector it will always return
2411: a device pointer to the device memory that contains this processor's portion of the vector data.
2413: Logically Collective; No Fortran Support
2415: Input Parameter:
2416: . x - the vector
2418: Output Parameters:
2419: + a - the array
2420: - mtype - memory type of the array
2422: Level: beginner
2424: Note:
2425: The array must be returned using a matching call to `VecRestoreArrayWriteAndMemType()`, where it will label the device memory as most recent.
2427: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2428: @*/
2429: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2430: {
2431: PetscFunctionBegin;
2434: PetscCall(VecSetErrorIfLocked(x, 1));
2435: PetscAssertPointer(a, 2);
2436: if (mtype) PetscAssertPointer(mtype, 3);
2437: if (x->ops->getarraywriteandmemtype) {
2438: /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2439: PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2440: } else if (x->ops->getarrayandmemtype) {
2441: PetscCall(VecGetArrayAndMemType(x, a, mtype));
2442: } else {
2443: /* VECNEST, VECVIENNACL */
2444: PetscCall(VecGetArrayWrite(x, a));
2445: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2446: }
2447: PetscFunctionReturn(PETSC_SUCCESS);
2448: }
2450: /*@C
2451: VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`
2453: Logically Collective; No Fortran Support
2455: Input Parameters:
2456: + x - the vector
2457: - a - the array
2459: Level: beginner
2461: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2462: @*/
2463: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar **a)
2464: {
2465: PetscFunctionBegin;
2468: PetscCall(VecSetErrorIfLocked(x, 1));
2469: if (a) PetscAssertPointer(a, 2);
2470: if (x->ops->restorearraywriteandmemtype) {
2471: /* VECCUDA/VECHIP */
2472: PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2473: PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2474: } else if (x->ops->restorearrayandmemtype) {
2475: PetscCall(VecRestoreArrayAndMemType(x, a));
2476: } else {
2477: PetscCall(VecRestoreArray(x, a));
2478: }
2479: if (a) *a = NULL;
2480: PetscFunctionReturn(PETSC_SUCCESS);
2481: }
2483: /*@
2484: VecPlaceArray - Allows one to replace the array in a vector with an
2485: array provided by the user. This is useful to avoid copying an array
2486: into a vector.
2488: Logically Collective; No Fortran Support
2490: Input Parameters:
2491: + vec - the vector
2492: - array - the array
2494: Level: developer
2496: Notes:
2497: Use `VecReplaceArray()` instead to permanently replace the array
2499: You can return to the original array with a call to `VecResetArray()`. `vec` does not take
2500: ownership of `array` in any way.
2502: The user must free `array` themselves but be careful not to
2503: do so before the vector has either been destroyed, had its original array restored with
2504: `VecResetArray()` or permanently replaced with `VecReplaceArray()`.
2506: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2507: @*/
2508: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2509: {
2510: PetscFunctionBegin;
2513: if (array) PetscAssertPointer(array, 2);
2514: PetscUseTypeMethod(vec, placearray, array);
2515: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2516: PetscFunctionReturn(PETSC_SUCCESS);
2517: }
2519: /*@C
2520: VecReplaceArray - Allows one to replace the array in a vector with an
2521: array provided by the user. This is useful to avoid copying an array
2522: into a vector.
2524: Logically Collective; No Fortran Support
2526: Input Parameters:
2527: + vec - the vector
2528: - array - the array
2530: Level: developer
2532: Notes:
2533: This permanently replaces the array and frees the memory associated
2534: with the old array. Use `VecPlaceArray()` to temporarily replace the array.
2536: The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
2537: freed by the user. It will be freed when the vector is destroyed.
2539: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2540: @*/
2541: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2542: {
2543: PetscFunctionBegin;
2546: PetscUseTypeMethod(vec, replacearray, array);
2547: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2548: PetscFunctionReturn(PETSC_SUCCESS);
2549: }
2551: /*MC
2552: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
2553: and makes them accessible via a Fortran pointer.
2555: Synopsis:
2556: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
2558: Collective
2560: Input Parameters:
2561: + x - a vector to mimic
2562: - n - the number of vectors to obtain
2564: Output Parameters:
2565: + y - Fortran pointer to the array of vectors
2566: - ierr - error code
2568: Example of Usage:
2569: .vb
2570: #include <petsc/finclude/petscvec.h>
2571: use petscvec
2573: Vec x
2574: Vec, pointer :: y(:)
2575: ....
2576: call VecDuplicateVecsF90(x,2,y,ierr)
2577: call VecSet(y(2),alpha,ierr)
2578: call VecSet(y(2),alpha,ierr)
2579: ....
2580: call VecDestroyVecsF90(2,y,ierr)
2581: .ve
2583: Level: beginner
2585: Note:
2586: Use `VecDestroyVecsF90()` to free the space.
2588: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecsF90()`, `VecDuplicateVecs()`
2589: M*/
2591: /*MC
2592: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2593: `VecGetArrayF90()`.
2595: Synopsis:
2596: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2598: Logically Collective
2600: Input Parameters:
2601: + x - vector
2602: - xx_v - the Fortran pointer to the array
2604: Output Parameter:
2605: . ierr - error code
2607: Example of Usage:
2608: .vb
2609: #include <petsc/finclude/petscvec.h>
2610: use petscvec
2612: PetscScalar, pointer :: xx_v(:)
2613: ....
2614: call VecGetArrayF90(x,xx_v,ierr)
2615: xx_v(3) = a
2616: call VecRestoreArrayF90(x,xx_v,ierr)
2617: .ve
2619: Level: beginner
2621: .seealso: [](ch_vectors), `Vec`, `VecGetArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrayReadF90()`
2622: M*/
2624: /*MC
2625: VecDestroyVecsF90 - Frees a block of vectors obtained with `VecDuplicateVecsF90()`.
2627: Synopsis:
2628: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
2630: Collective
2632: Input Parameters:
2633: + n - the number of vectors previously obtained
2634: - x - pointer to array of vector pointers
2636: Output Parameter:
2637: . ierr - error code
2639: Level: beginner
2641: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecs()`, `VecDuplicateVecsF90()`
2642: M*/
2644: /*MC
2645: VecGetArrayF90 - Accesses a vector array from Fortran. For default PETSc
2646: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2647: this routine is implementation dependent. You MUST call `VecRestoreArrayF90()`
2648: when you no longer need access to the array.
2650: Synopsis:
2651: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2653: Logically Collective
2655: Input Parameter:
2656: . x - vector
2658: Output Parameters:
2659: + xx_v - the Fortran pointer to the array
2660: - ierr - error code
2662: Example of Usage:
2663: .vb
2664: #include <petsc/finclude/petscvec.h>
2665: use petscvec
2667: PetscScalar, pointer :: xx_v(:)
2668: ....
2669: call VecGetArrayF90(x,xx_v,ierr)
2670: xx_v(3) = a
2671: call VecRestoreArrayF90(x,xx_v,ierr)
2672: .ve
2674: Level: beginner
2676: Note:
2677: If you ONLY intend to read entries from the array and not change any entries you should use `VecGetArrayReadF90()`.
2679: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayReadF90()`
2680: M*/
2682: /*MC
2683: VecGetArrayReadF90 - Accesses a read only array from Fortran. For default PETSc
2684: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2685: this routine is implementation dependent. You MUST call `VecRestoreArrayReadF90()`
2686: when you no longer need access to the array.
2688: Synopsis:
2689: VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2691: Logically Collective
2693: Input Parameter:
2694: . x - vector
2696: Output Parameters:
2697: + xx_v - the Fortran pointer to the array
2698: - ierr - error code
2700: Example of Usage:
2701: .vb
2702: #include <petsc/finclude/petscvec.h>
2703: use petscvec
2705: PetscScalar, pointer :: xx_v(:)
2706: ....
2707: call VecGetArrayReadF90(x,xx_v,ierr)
2708: a = xx_v(3)
2709: call VecRestoreArrayReadF90(x,xx_v,ierr)
2710: .ve
2712: Level: beginner
2714: Note:
2715: If you intend to write entries into the array you must use `VecGetArrayF90()`.
2717: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecGetArrayF90()`
2718: M*/
2720: /*MC
2721: VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2722: `VecGetArrayReadF90()`.
2724: Synopsis:
2725: VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2727: Logically Collective
2729: Input Parameters:
2730: + x - vector
2731: - xx_v - the Fortran pointer to the array
2733: Output Parameter:
2734: . ierr - error code
2736: Example of Usage:
2737: .vb
2738: #include <petsc/finclude/petscvec.h>
2739: use petscvec
2741: PetscScalar, pointer :: xx_v(:)
2742: ....
2743: call VecGetArrayReadF90(x,xx_v,ierr)
2744: a = xx_v(3)
2745: call VecRestoreArrayReadF90(x,xx_v,ierr)
2746: .ve
2748: Level: beginner
2750: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecRestoreArrayF90()`
2751: M*/
2753: /*@C
2754: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2755: processor's portion of the vector data. You MUST call `VecRestoreArray2d()`
2756: when you no longer need access to the array.
2758: Logically Collective
2760: Input Parameters:
2761: + x - the vector
2762: . m - first dimension of two dimensional array
2763: . n - second dimension of two dimensional array
2764: . mstart - first index you will use in first coordinate direction (often 0)
2765: - nstart - first index in the second coordinate direction (often 0)
2767: Output Parameter:
2768: . a - location to put pointer to the array
2770: Level: developer
2772: Notes:
2773: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2774: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2775: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2776: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2778: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2780: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2781: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2782: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2783: @*/
2784: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2785: {
2786: PetscInt i, N;
2787: PetscScalar *aa;
2789: PetscFunctionBegin;
2791: PetscAssertPointer(a, 6);
2793: PetscCall(VecGetLocalSize(x, &N));
2794: 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);
2795: PetscCall(VecGetArray(x, &aa));
2797: PetscCall(PetscMalloc1(m, a));
2798: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2799: *a -= mstart;
2800: PetscFunctionReturn(PETSC_SUCCESS);
2801: }
2803: /*@C
2804: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2805: processor's portion of the vector data. You MUST call `VecRestoreArray2dWrite()`
2806: when you no longer need access to the array.
2808: Logically Collective
2810: Input Parameters:
2811: + x - the vector
2812: . m - first dimension of two dimensional array
2813: . n - second dimension of two dimensional array
2814: . mstart - first index you will use in first coordinate direction (often 0)
2815: - nstart - first index in the second coordinate direction (often 0)
2817: Output Parameter:
2818: . a - location to put pointer to the array
2820: Level: developer
2822: Notes:
2823: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2824: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2825: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2826: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2828: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2830: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2831: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2832: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2833: @*/
2834: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2835: {
2836: PetscInt i, N;
2837: PetscScalar *aa;
2839: PetscFunctionBegin;
2841: PetscAssertPointer(a, 6);
2843: PetscCall(VecGetLocalSize(x, &N));
2844: 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);
2845: PetscCall(VecGetArrayWrite(x, &aa));
2847: PetscCall(PetscMalloc1(m, a));
2848: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2849: *a -= mstart;
2850: PetscFunctionReturn(PETSC_SUCCESS);
2851: }
2853: /*@C
2854: VecRestoreArray2d - Restores a vector after `VecGetArray2d()` has been called.
2856: Logically Collective
2858: Input Parameters:
2859: + x - the vector
2860: . m - first dimension of two dimensional array
2861: . n - second dimension of the two dimensional array
2862: . mstart - first index you will use in first coordinate direction (often 0)
2863: . nstart - first index in the second coordinate direction (often 0)
2864: - a - location of pointer to array obtained from `VecGetArray2d()`
2866: Level: developer
2868: Notes:
2869: For regular PETSc vectors this routine does not involve any copies. For
2870: any special vectors that do not store local vector data in a contiguous
2871: array, this routine will copy the data back into the underlying
2872: vector data structure from the array obtained with `VecGetArray()`.
2874: This routine actually zeros out the `a` pointer.
2876: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2877: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2878: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2879: @*/
2880: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2881: {
2882: void *dummy;
2884: PetscFunctionBegin;
2886: PetscAssertPointer(a, 6);
2888: dummy = (void *)(*a + mstart);
2889: PetscCall(PetscFree(dummy));
2890: PetscCall(VecRestoreArray(x, NULL));
2891: *a = NULL;
2892: PetscFunctionReturn(PETSC_SUCCESS);
2893: }
2895: /*@C
2896: VecRestoreArray2dWrite - Restores a vector after `VecGetArray2dWrite()` has been called.
2898: Logically Collective
2900: Input Parameters:
2901: + x - the vector
2902: . m - first dimension of two dimensional array
2903: . n - second dimension of the two dimensional array
2904: . mstart - first index you will use in first coordinate direction (often 0)
2905: . nstart - first index in the second coordinate direction (often 0)
2906: - a - location of pointer to array obtained from `VecGetArray2d()`
2908: Level: developer
2910: Notes:
2911: For regular PETSc vectors this routine does not involve any copies. For
2912: any special vectors that do not store local vector data in a contiguous
2913: array, this routine will copy the data back into the underlying
2914: vector data structure from the array obtained with `VecGetArray()`.
2916: This routine actually zeros out the `a` pointer.
2918: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2919: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2920: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2921: @*/
2922: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2923: {
2924: void *dummy;
2926: PetscFunctionBegin;
2928: PetscAssertPointer(a, 6);
2930: dummy = (void *)(*a + mstart);
2931: PetscCall(PetscFree(dummy));
2932: PetscCall(VecRestoreArrayWrite(x, NULL));
2933: PetscFunctionReturn(PETSC_SUCCESS);
2934: }
2936: /*@C
2937: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2938: processor's portion of the vector data. You MUST call `VecRestoreArray1d()`
2939: when you no longer need access to the array.
2941: Logically Collective
2943: Input Parameters:
2944: + x - the vector
2945: . m - first dimension of two dimensional array
2946: - mstart - first index you will use in first coordinate direction (often 0)
2948: Output Parameter:
2949: . a - location to put pointer to the array
2951: Level: developer
2953: Notes:
2954: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2955: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2956: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2958: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2960: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2961: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2962: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2963: @*/
2964: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2965: {
2966: PetscInt N;
2968: PetscFunctionBegin;
2970: PetscAssertPointer(a, 4);
2972: PetscCall(VecGetLocalSize(x, &N));
2973: 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);
2974: PetscCall(VecGetArray(x, a));
2975: *a -= mstart;
2976: PetscFunctionReturn(PETSC_SUCCESS);
2977: }
2979: /*@C
2980: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
2981: processor's portion of the vector data. You MUST call `VecRestoreArray1dWrite()`
2982: when you no longer need access to the array.
2984: Logically Collective
2986: Input Parameters:
2987: + x - the vector
2988: . m - first dimension of two dimensional array
2989: - mstart - first index you will use in first coordinate direction (often 0)
2991: Output Parameter:
2992: . a - location to put pointer to the array
2994: Level: developer
2996: Notes:
2997: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2998: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2999: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3001: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3003: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3004: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3005: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3006: @*/
3007: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3008: {
3009: PetscInt N;
3011: PetscFunctionBegin;
3013: PetscAssertPointer(a, 4);
3015: PetscCall(VecGetLocalSize(x, &N));
3016: 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);
3017: PetscCall(VecGetArrayWrite(x, a));
3018: *a -= mstart;
3019: PetscFunctionReturn(PETSC_SUCCESS);
3020: }
3022: /*@C
3023: VecRestoreArray1d - Restores a vector after `VecGetArray1d()` has been called.
3025: Logically Collective
3027: Input Parameters:
3028: + x - the vector
3029: . m - first dimension of two dimensional array
3030: . mstart - first index you will use in first coordinate direction (often 0)
3031: - a - location of pointer to array obtained from `VecGetArray1d()`
3033: Level: developer
3035: Notes:
3036: For regular PETSc vectors this routine does not involve any copies. For
3037: any special vectors that do not store local vector data in a contiguous
3038: array, this routine will copy the data back into the underlying
3039: vector data structure from the array obtained with `VecGetArray1d()`.
3041: This routine actually zeros out the `a` pointer.
3043: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3044: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3045: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3046: @*/
3047: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3048: {
3049: PetscFunctionBegin;
3052: PetscCall(VecRestoreArray(x, NULL));
3053: *a = NULL;
3054: PetscFunctionReturn(PETSC_SUCCESS);
3055: }
3057: /*@C
3058: VecRestoreArray1dWrite - Restores a vector after `VecGetArray1dWrite()` has been called.
3060: Logically Collective
3062: Input Parameters:
3063: + x - the vector
3064: . m - first dimension of two dimensional array
3065: . mstart - first index you will use in first coordinate direction (often 0)
3066: - a - location of pointer to array obtained from `VecGetArray1d()`
3068: Level: developer
3070: Notes:
3071: For regular PETSc vectors this routine does not involve any copies. For
3072: any special vectors that do not store local vector data in a contiguous
3073: array, this routine will copy the data back into the underlying
3074: vector data structure from the array obtained with `VecGetArray1d()`.
3076: This routine actually zeros out the `a` pointer.
3078: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3079: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3080: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3081: @*/
3082: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3083: {
3084: PetscFunctionBegin;
3087: PetscCall(VecRestoreArrayWrite(x, NULL));
3088: *a = NULL;
3089: PetscFunctionReturn(PETSC_SUCCESS);
3090: }
3092: /*@C
3093: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
3094: processor's portion of the vector data. You MUST call `VecRestoreArray3d()`
3095: when you no longer need access to the array.
3097: Logically Collective
3099: Input Parameters:
3100: + x - the vector
3101: . m - first dimension of three dimensional array
3102: . n - second dimension of three dimensional array
3103: . p - third dimension of three dimensional array
3104: . mstart - first index you will use in first coordinate direction (often 0)
3105: . nstart - first index in the second coordinate direction (often 0)
3106: - pstart - first index in the third coordinate direction (often 0)
3108: Output Parameter:
3109: . a - location to put pointer to the array
3111: Level: developer
3113: Notes:
3114: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3115: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3116: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3117: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3119: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3121: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3122: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
3123: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3124: @*/
3125: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3126: {
3127: PetscInt i, N, j;
3128: PetscScalar *aa, **b;
3130: PetscFunctionBegin;
3132: PetscAssertPointer(a, 8);
3134: PetscCall(VecGetLocalSize(x, &N));
3135: 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);
3136: PetscCall(VecGetArray(x, &aa));
3138: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3139: b = (PetscScalar **)((*a) + m);
3140: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3141: for (i = 0; i < m; i++)
3142: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3143: *a -= mstart;
3144: PetscFunctionReturn(PETSC_SUCCESS);
3145: }
3147: /*@C
3148: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
3149: processor's portion of the vector data. You MUST call `VecRestoreArray3dWrite()`
3150: when you no longer need access to the array.
3152: Logically Collective
3154: Input Parameters:
3155: + x - the vector
3156: . m - first dimension of three dimensional array
3157: . n - second dimension of three dimensional array
3158: . p - third dimension of three dimensional array
3159: . mstart - first index you will use in first coordinate direction (often 0)
3160: . nstart - first index in the second coordinate direction (often 0)
3161: - pstart - first index in the third coordinate direction (often 0)
3163: Output Parameter:
3164: . a - location to put pointer to the array
3166: Level: developer
3168: Notes:
3169: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3170: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3171: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3172: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3174: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3176: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3177: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3178: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3179: @*/
3180: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3181: {
3182: PetscInt i, N, j;
3183: PetscScalar *aa, **b;
3185: PetscFunctionBegin;
3187: PetscAssertPointer(a, 8);
3189: PetscCall(VecGetLocalSize(x, &N));
3190: 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);
3191: PetscCall(VecGetArrayWrite(x, &aa));
3193: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3194: b = (PetscScalar **)((*a) + m);
3195: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3196: for (i = 0; i < m; i++)
3197: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3199: *a -= mstart;
3200: PetscFunctionReturn(PETSC_SUCCESS);
3201: }
3203: /*@C
3204: VecRestoreArray3d - Restores a vector after `VecGetArray3d()` has been called.
3206: Logically Collective
3208: Input Parameters:
3209: + x - the vector
3210: . m - first dimension of three dimensional array
3211: . n - second dimension of the three dimensional array
3212: . p - third dimension of the three dimensional array
3213: . mstart - first index you will use in first coordinate direction (often 0)
3214: . nstart - first index in the second coordinate direction (often 0)
3215: . pstart - first index in the third coordinate direction (often 0)
3216: - a - location of pointer to array obtained from VecGetArray3d()
3218: Level: developer
3220: Notes:
3221: For regular PETSc vectors this routine does not involve any copies. For
3222: any special vectors that do not store local vector data in a contiguous
3223: array, this routine will copy the data back into the underlying
3224: vector data structure from the array obtained with `VecGetArray()`.
3226: This routine actually zeros out the `a` pointer.
3228: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3229: `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3230: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3231: @*/
3232: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3233: {
3234: void *dummy;
3236: PetscFunctionBegin;
3238: PetscAssertPointer(a, 8);
3240: dummy = (void *)(*a + mstart);
3241: PetscCall(PetscFree(dummy));
3242: PetscCall(VecRestoreArray(x, NULL));
3243: *a = NULL;
3244: PetscFunctionReturn(PETSC_SUCCESS);
3245: }
3247: /*@C
3248: VecRestoreArray3dWrite - Restores a vector after `VecGetArray3dWrite()` has been called.
3250: Logically Collective
3252: Input Parameters:
3253: + x - the vector
3254: . m - first dimension of three dimensional array
3255: . n - second dimension of the three dimensional array
3256: . p - third dimension of the three dimensional array
3257: . mstart - first index you will use in first coordinate direction (often 0)
3258: . nstart - first index in the second coordinate direction (often 0)
3259: . pstart - first index in the third coordinate direction (often 0)
3260: - a - location of pointer to array obtained from VecGetArray3d()
3262: Level: developer
3264: Notes:
3265: For regular PETSc vectors this routine does not involve any copies. For
3266: any special vectors that do not store local vector data in a contiguous
3267: array, this routine will copy the data back into the underlying
3268: vector data structure from the array obtained with `VecGetArray()`.
3270: This routine actually zeros out the `a` pointer.
3272: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3273: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3274: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3275: @*/
3276: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3277: {
3278: void *dummy;
3280: PetscFunctionBegin;
3282: PetscAssertPointer(a, 8);
3284: dummy = (void *)(*a + mstart);
3285: PetscCall(PetscFree(dummy));
3286: PetscCall(VecRestoreArrayWrite(x, NULL));
3287: *a = NULL;
3288: PetscFunctionReturn(PETSC_SUCCESS);
3289: }
3291: /*@C
3292: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3293: processor's portion of the vector data. You MUST call `VecRestoreArray4d()`
3294: when you no longer need access to the array.
3296: Logically Collective
3298: Input Parameters:
3299: + x - the vector
3300: . m - first dimension of four dimensional array
3301: . n - second dimension of four dimensional array
3302: . p - third dimension of four dimensional array
3303: . q - fourth dimension of four dimensional array
3304: . mstart - first index you will use in first coordinate direction (often 0)
3305: . nstart - first index in the second coordinate direction (often 0)
3306: . pstart - first index in the third coordinate direction (often 0)
3307: - qstart - first index in the fourth coordinate direction (often 0)
3309: Output Parameter:
3310: . a - location to put pointer to the array
3312: Level: developer
3314: Notes:
3315: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3316: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3317: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3318: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3320: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3322: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3323: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3324: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3325: @*/
3326: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3327: {
3328: PetscInt i, N, j, k;
3329: PetscScalar *aa, ***b, **c;
3331: PetscFunctionBegin;
3333: PetscAssertPointer(a, 10);
3335: PetscCall(VecGetLocalSize(x, &N));
3336: 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);
3337: PetscCall(VecGetArray(x, &aa));
3339: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3340: b = (PetscScalar ***)((*a) + m);
3341: c = (PetscScalar **)(b + m * n);
3342: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3343: for (i = 0; i < m; i++)
3344: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3345: for (i = 0; i < m; i++)
3346: for (j = 0; j < n; j++)
3347: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3348: *a -= mstart;
3349: PetscFunctionReturn(PETSC_SUCCESS);
3350: }
3352: /*@C
3353: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3354: processor's portion of the vector data. You MUST call `VecRestoreArray4dWrite()`
3355: when you no longer need access to the array.
3357: Logically Collective
3359: Input Parameters:
3360: + x - the vector
3361: . m - first dimension of four dimensional array
3362: . n - second dimension of four dimensional array
3363: . p - third dimension of four dimensional array
3364: . q - fourth dimension of four dimensional array
3365: . mstart - first index you will use in first coordinate direction (often 0)
3366: . nstart - first index in the second coordinate direction (often 0)
3367: . pstart - first index in the third coordinate direction (often 0)
3368: - qstart - first index in the fourth coordinate direction (often 0)
3370: Output Parameter:
3371: . a - location to put pointer to the array
3373: Level: developer
3375: Notes:
3376: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3377: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3378: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3379: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3381: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3383: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3384: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3385: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3386: @*/
3387: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3388: {
3389: PetscInt i, N, j, k;
3390: PetscScalar *aa, ***b, **c;
3392: PetscFunctionBegin;
3394: PetscAssertPointer(a, 10);
3396: PetscCall(VecGetLocalSize(x, &N));
3397: 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);
3398: PetscCall(VecGetArrayWrite(x, &aa));
3400: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3401: b = (PetscScalar ***)((*a) + m);
3402: c = (PetscScalar **)(b + m * n);
3403: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3404: for (i = 0; i < m; i++)
3405: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3406: for (i = 0; i < m; i++)
3407: for (j = 0; j < n; j++)
3408: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3409: *a -= mstart;
3410: PetscFunctionReturn(PETSC_SUCCESS);
3411: }
3413: /*@C
3414: VecRestoreArray4d - Restores a vector after `VecGetArray4d()` has been called.
3416: Logically Collective
3418: Input Parameters:
3419: + x - the vector
3420: . m - first dimension of four dimensional array
3421: . n - second dimension of the four dimensional array
3422: . p - third dimension of the four dimensional array
3423: . q - fourth dimension of the four dimensional array
3424: . mstart - first index you will use in first coordinate direction (often 0)
3425: . nstart - first index in the second coordinate direction (often 0)
3426: . pstart - first index in the third coordinate direction (often 0)
3427: . qstart - first index in the fourth coordinate direction (often 0)
3428: - a - location of pointer to array obtained from VecGetArray4d()
3430: Level: developer
3432: Notes:
3433: For regular PETSc vectors this routine does not involve any copies. For
3434: any special vectors that do not store local vector data in a contiguous
3435: array, this routine will copy the data back into the underlying
3436: vector data structure from the array obtained with `VecGetArray()`.
3438: This routine actually zeros out the `a` pointer.
3440: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3441: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3442: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecGet`
3443: @*/
3444: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3445: {
3446: void *dummy;
3448: PetscFunctionBegin;
3450: PetscAssertPointer(a, 10);
3452: dummy = (void *)(*a + mstart);
3453: PetscCall(PetscFree(dummy));
3454: PetscCall(VecRestoreArray(x, NULL));
3455: *a = NULL;
3456: PetscFunctionReturn(PETSC_SUCCESS);
3457: }
3459: /*@C
3460: VecRestoreArray4dWrite - Restores a vector after `VecGetArray4dWrite()` has been called.
3462: Logically Collective
3464: Input Parameters:
3465: + x - the vector
3466: . m - first dimension of four dimensional array
3467: . n - second dimension of the four dimensional array
3468: . p - third dimension of the four dimensional array
3469: . q - fourth dimension of the four dimensional array
3470: . mstart - first index you will use in first coordinate direction (often 0)
3471: . nstart - first index in the second coordinate direction (often 0)
3472: . pstart - first index in the third coordinate direction (often 0)
3473: . qstart - first index in the fourth coordinate direction (often 0)
3474: - a - location of pointer to array obtained from `VecGetArray4d()`
3476: Level: developer
3478: Notes:
3479: For regular PETSc vectors this routine does not involve any copies. For
3480: any special vectors that do not store local vector data in a contiguous
3481: array, this routine will copy the data back into the underlying
3482: vector data structure from the array obtained with `VecGetArray()`.
3484: This routine actually zeros out the `a` pointer.
3486: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3487: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3488: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3489: @*/
3490: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3491: {
3492: void *dummy;
3494: PetscFunctionBegin;
3496: PetscAssertPointer(a, 10);
3498: dummy = (void *)(*a + mstart);
3499: PetscCall(PetscFree(dummy));
3500: PetscCall(VecRestoreArrayWrite(x, NULL));
3501: *a = NULL;
3502: PetscFunctionReturn(PETSC_SUCCESS);
3503: }
3505: /*@C
3506: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3507: processor's portion of the vector data. You MUST call `VecRestoreArray2dRead()`
3508: when you no longer need access to the array.
3510: Logically Collective
3512: Input Parameters:
3513: + x - the vector
3514: . m - first dimension of two dimensional array
3515: . n - second dimension of two dimensional array
3516: . mstart - first index you will use in first coordinate direction (often 0)
3517: - nstart - first index in the second coordinate direction (often 0)
3519: Output Parameter:
3520: . a - location to put pointer to the array
3522: Level: developer
3524: Notes:
3525: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
3526: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3527: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3528: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
3530: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3532: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3533: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3534: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3535: @*/
3536: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3537: {
3538: PetscInt i, N;
3539: const PetscScalar *aa;
3541: PetscFunctionBegin;
3543: PetscAssertPointer(a, 6);
3545: PetscCall(VecGetLocalSize(x, &N));
3546: 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);
3547: PetscCall(VecGetArrayRead(x, &aa));
3549: PetscCall(PetscMalloc1(m, a));
3550: for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3551: *a -= mstart;
3552: PetscFunctionReturn(PETSC_SUCCESS);
3553: }
3555: /*@C
3556: VecRestoreArray2dRead - Restores a vector after `VecGetArray2dRead()` has been called.
3558: Logically Collective
3560: Input Parameters:
3561: + x - the vector
3562: . m - first dimension of two dimensional array
3563: . n - second dimension of the two dimensional array
3564: . mstart - first index you will use in first coordinate direction (often 0)
3565: . nstart - first index in the second coordinate direction (often 0)
3566: - a - location of pointer to array obtained from VecGetArray2d()
3568: Level: developer
3570: Notes:
3571: For regular PETSc vectors this routine does not involve any copies. For
3572: any special vectors that do not store local vector data in a contiguous
3573: array, this routine will copy the data back into the underlying
3574: vector data structure from the array obtained with `VecGetArray()`.
3576: This routine actually zeros out the `a` pointer.
3578: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3579: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3580: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3581: @*/
3582: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3583: {
3584: void *dummy;
3586: PetscFunctionBegin;
3588: PetscAssertPointer(a, 6);
3590: dummy = (void *)(*a + mstart);
3591: PetscCall(PetscFree(dummy));
3592: PetscCall(VecRestoreArrayRead(x, NULL));
3593: *a = NULL;
3594: PetscFunctionReturn(PETSC_SUCCESS);
3595: }
3597: /*@C
3598: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3599: processor's portion of the vector data. You MUST call `VecRestoreArray1dRead()`
3600: when you no longer need access to the array.
3602: Logically Collective
3604: Input Parameters:
3605: + x - the vector
3606: . m - first dimension of two dimensional array
3607: - mstart - first index you will use in first coordinate direction (often 0)
3609: Output Parameter:
3610: . a - location to put pointer to the array
3612: Level: developer
3614: Notes:
3615: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3616: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3617: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3619: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3621: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3622: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3623: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3624: @*/
3625: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3626: {
3627: PetscInt N;
3629: PetscFunctionBegin;
3631: PetscAssertPointer(a, 4);
3633: PetscCall(VecGetLocalSize(x, &N));
3634: 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);
3635: PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3636: *a -= mstart;
3637: PetscFunctionReturn(PETSC_SUCCESS);
3638: }
3640: /*@C
3641: VecRestoreArray1dRead - Restores a vector after `VecGetArray1dRead()` has been called.
3643: Logically Collective
3645: Input Parameters:
3646: + x - the vector
3647: . m - first dimension of two dimensional array
3648: . mstart - first index you will use in first coordinate direction (often 0)
3649: - a - location of pointer to array obtained from `VecGetArray1dRead()`
3651: Level: developer
3653: Notes:
3654: For regular PETSc vectors this routine does not involve any copies. For
3655: any special vectors that do not store local vector data in a contiguous
3656: array, this routine will copy the data back into the underlying
3657: vector data structure from the array obtained with `VecGetArray1dRead()`.
3659: This routine actually zeros out the `a` pointer.
3661: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3662: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3663: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3664: @*/
3665: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3666: {
3667: PetscFunctionBegin;
3670: PetscCall(VecRestoreArrayRead(x, NULL));
3671: *a = NULL;
3672: PetscFunctionReturn(PETSC_SUCCESS);
3673: }
3675: /*@C
3676: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3677: processor's portion of the vector data. You MUST call `VecRestoreArray3dRead()`
3678: when you no longer need access to the array.
3680: Logically Collective
3682: Input Parameters:
3683: + x - the vector
3684: . m - first dimension of three dimensional array
3685: . n - second dimension of three dimensional array
3686: . p - third dimension of three dimensional array
3687: . mstart - first index you will use in first coordinate direction (often 0)
3688: . nstart - first index in the second coordinate direction (often 0)
3689: - pstart - first index in the third coordinate direction (often 0)
3691: Output Parameter:
3692: . a - location to put pointer to the array
3694: Level: developer
3696: Notes:
3697: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3698: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3699: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3700: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3dRead()`.
3702: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3704: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3705: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3706: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3707: @*/
3708: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3709: {
3710: PetscInt i, N, j;
3711: const PetscScalar *aa;
3712: PetscScalar **b;
3714: PetscFunctionBegin;
3716: PetscAssertPointer(a, 8);
3718: PetscCall(VecGetLocalSize(x, &N));
3719: 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);
3720: PetscCall(VecGetArrayRead(x, &aa));
3722: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3723: b = (PetscScalar **)((*a) + m);
3724: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3725: for (i = 0; i < m; i++)
3726: for (j = 0; j < n; j++) b[i * n + j] = (PetscScalar *)aa + i * n * p + j * p - pstart;
3727: *a -= mstart;
3728: PetscFunctionReturn(PETSC_SUCCESS);
3729: }
3731: /*@C
3732: VecRestoreArray3dRead - Restores a vector after `VecGetArray3dRead()` has been called.
3734: Logically Collective
3736: Input Parameters:
3737: + x - the vector
3738: . m - first dimension of three dimensional array
3739: . n - second dimension of the three dimensional array
3740: . p - third dimension of the three dimensional array
3741: . mstart - first index you will use in first coordinate direction (often 0)
3742: . nstart - first index in the second coordinate direction (often 0)
3743: . pstart - first index in the third coordinate direction (often 0)
3744: - a - location of pointer to array obtained from `VecGetArray3dRead()`
3746: Level: developer
3748: Notes:
3749: For regular PETSc vectors this routine does not involve any copies. For
3750: any special vectors that do not store local vector data in a contiguous
3751: array, this routine will copy the data back into the underlying
3752: vector data structure from the array obtained with `VecGetArray()`.
3754: This routine actually zeros out the `a` pointer.
3756: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3757: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3758: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3759: @*/
3760: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3761: {
3762: void *dummy;
3764: PetscFunctionBegin;
3766: PetscAssertPointer(a, 8);
3768: dummy = (void *)(*a + mstart);
3769: PetscCall(PetscFree(dummy));
3770: PetscCall(VecRestoreArrayRead(x, NULL));
3771: *a = NULL;
3772: PetscFunctionReturn(PETSC_SUCCESS);
3773: }
3775: /*@C
3776: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3777: processor's portion of the vector data. You MUST call `VecRestoreArray4dRead()`
3778: when you no longer need access to the array.
3780: Logically Collective
3782: Input Parameters:
3783: + x - the vector
3784: . m - first dimension of four dimensional array
3785: . n - second dimension of four dimensional array
3786: . p - third dimension of four dimensional array
3787: . q - fourth dimension of four dimensional array
3788: . mstart - first index you will use in first coordinate direction (often 0)
3789: . nstart - first index in the second coordinate direction (often 0)
3790: . pstart - first index in the third coordinate direction (often 0)
3791: - qstart - first index in the fourth coordinate direction (often 0)
3793: Output Parameter:
3794: . a - location to put pointer to the array
3796: Level: beginner
3798: Notes:
3799: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3800: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3801: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3802: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3804: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3806: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3807: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3808: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3809: @*/
3810: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3811: {
3812: PetscInt i, N, j, k;
3813: const PetscScalar *aa;
3814: PetscScalar ***b, **c;
3816: PetscFunctionBegin;
3818: PetscAssertPointer(a, 10);
3820: PetscCall(VecGetLocalSize(x, &N));
3821: 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);
3822: PetscCall(VecGetArrayRead(x, &aa));
3824: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3825: b = (PetscScalar ***)((*a) + m);
3826: c = (PetscScalar **)(b + m * n);
3827: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3828: for (i = 0; i < m; i++)
3829: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3830: for (i = 0; i < m; i++)
3831: for (j = 0; j < n; j++)
3832: 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;
3833: *a -= mstart;
3834: PetscFunctionReturn(PETSC_SUCCESS);
3835: }
3837: /*@C
3838: VecRestoreArray4dRead - Restores a vector after `VecGetArray4d()` has been called.
3840: Logically Collective
3842: Input Parameters:
3843: + x - the vector
3844: . m - first dimension of four dimensional array
3845: . n - second dimension of the four dimensional array
3846: . p - third dimension of the four dimensional array
3847: . q - fourth dimension of the four dimensional array
3848: . mstart - first index you will use in first coordinate direction (often 0)
3849: . nstart - first index in the second coordinate direction (often 0)
3850: . pstart - first index in the third coordinate direction (often 0)
3851: . qstart - first index in the fourth coordinate direction (often 0)
3852: - a - location of pointer to array obtained from `VecGetArray4dRead()`
3854: Level: beginner
3856: Notes:
3857: For regular PETSc vectors this routine does not involve any copies. For
3858: any special vectors that do not store local vector data in a contiguous
3859: array, this routine will copy the data back into the underlying
3860: vector data structure from the array obtained with `VecGetArray()`.
3862: This routine actually zeros out the `a` pointer.
3864: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3865: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3866: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3867: @*/
3868: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3869: {
3870: void *dummy;
3872: PetscFunctionBegin;
3874: PetscAssertPointer(a, 10);
3876: dummy = (void *)(*a + mstart);
3877: PetscCall(PetscFree(dummy));
3878: PetscCall(VecRestoreArrayRead(x, NULL));
3879: *a = NULL;
3880: PetscFunctionReturn(PETSC_SUCCESS);
3881: }
3883: #if defined(PETSC_USE_DEBUG)
3885: /*@
3886: VecLockGet - Gets the current lock status of a vector
3888: Logically Collective
3890: Input Parameter:
3891: . x - the vector
3893: Output Parameter:
3894: . state - greater than zero indicates the vector is locked for read; less than zero indicates the vector is
3895: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3897: Level: advanced
3899: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3900: @*/
3901: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3902: {
3903: PetscFunctionBegin;
3905: *state = x->lock;
3906: PetscFunctionReturn(PETSC_SUCCESS);
3907: }
3909: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3910: {
3911: PetscFunctionBegin;
3913: PetscAssertPointer(file, 2);
3914: PetscAssertPointer(func, 3);
3915: PetscAssertPointer(line, 4);
3916: #if !PetscDefined(HAVE_THREADSAFETY)
3917: {
3918: const int index = x->lockstack.currentsize - 1;
3920: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupted vec lock stack, have negative index %d", index);
3921: *file = x->lockstack.file[index];
3922: *func = x->lockstack.function[index];
3923: *line = x->lockstack.line[index];
3924: }
3925: #else
3926: *file = NULL;
3927: *func = NULL;
3928: *line = 0;
3929: #endif
3930: PetscFunctionReturn(PETSC_SUCCESS);
3931: }
3933: /*@
3934: VecLockReadPush - Pushes a read-only lock on a vector to prevent it from being written to
3936: Logically Collective
3938: Input Parameter:
3939: . x - the vector
3941: Level: intermediate
3943: Notes:
3944: If this is set then calls to `VecGetArray()` or `VecSetValues()` or any other routines that change the vectors values will generate an error.
3946: The call can be nested, i.e., called multiple times on the same vector, but each `VecLockReadPush()` has to have one matching
3947: `VecLockReadPop()`, which removes the latest read-only lock.
3949: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3950: @*/
3951: PetscErrorCode VecLockReadPush(Vec x)
3952: {
3953: PetscFunctionBegin;
3955: 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");
3956: #if !PetscDefined(HAVE_THREADSAFETY)
3957: {
3958: const char *file, *func;
3959: int index, line;
3961: if ((index = petscstack.currentsize - 2) == -1) {
3962: // vec was locked "outside" of petsc, either in user-land or main. the error message will
3963: // now show this function as the culprit, but it will include the stacktrace
3964: file = "unknown user-file";
3965: func = "unknown_user_function";
3966: line = 0;
3967: } else {
3968: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected petscstack, have negative index %d", index);
3969: file = petscstack.file[index];
3970: func = petscstack.function[index];
3971: line = petscstack.line[index];
3972: }
3973: PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
3974: }
3975: #endif
3976: PetscFunctionReturn(PETSC_SUCCESS);
3977: }
3979: /*@
3980: VecLockReadPop - Pops a read-only lock from a vector
3982: Logically Collective
3984: Input Parameter:
3985: . x - the vector
3987: Level: intermediate
3989: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
3990: @*/
3991: PetscErrorCode VecLockReadPop(Vec x)
3992: {
3993: PetscFunctionBegin;
3995: PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
3996: #if !PetscDefined(HAVE_THREADSAFETY)
3997: {
3998: const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];
4000: PetscStackPop_Private(x->lockstack, previous);
4001: }
4002: #endif
4003: PetscFunctionReturn(PETSC_SUCCESS);
4004: }
4006: /*@C
4007: VecLockWriteSet - Lock or unlock a vector for exclusive read/write access
4009: Logically Collective
4011: Input Parameters:
4012: + x - the vector
4013: - flg - `PETSC_TRUE` to lock the vector for exclusive read/write access; `PETSC_FALSE` to unlock it.
4015: Level: intermediate
4017: Notes:
4018: The function is useful in split-phase computations, which usually have a begin phase and an end phase.
4019: One can call `VecLockWriteSet`(x,`PETSC_TRUE`) in the begin phase to lock a vector for exclusive
4020: access, and call `VecLockWriteSet`(x,`PETSC_FALSE`) in the end phase to unlock the vector from exclusive
4021: access. In this way, one is ensured no other operations can access the vector in between. The code may like
4023: .vb
4024: VecGetArray(x,&xdata); // begin phase
4025: VecLockWriteSet(v,PETSC_TRUE);
4027: Other operations, which can not access x anymore (they can access xdata, of course)
4029: VecRestoreArray(x,&vdata); // end phase
4030: VecLockWriteSet(v,PETSC_FALSE);
4031: .ve
4033: The call can not be nested on the same vector, in other words, one can not call `VecLockWriteSet`(x,`PETSC_TRUE`)
4034: again before calling `VecLockWriteSet`(v,`PETSC_FALSE`).
4036: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
4037: @*/
4038: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
4039: {
4040: PetscFunctionBegin;
4042: if (flg) {
4043: 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");
4044: 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");
4045: x->lock = -1;
4046: } else {
4047: 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");
4048: x->lock = 0;
4049: }
4050: PetscFunctionReturn(PETSC_SUCCESS);
4051: }
4053: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
4054: /*@
4055: VecLockPush - Pushes a read-only lock on a vector to prevent it from being written to
4057: Level: deprecated
4059: .seealso: [](ch_vectors), `Vec`, `VecLockReadPush()`
4060: @*/
4061: PetscErrorCode VecLockPush(Vec x)
4062: {
4063: PetscFunctionBegin;
4064: PetscCall(VecLockReadPush(x));
4065: PetscFunctionReturn(PETSC_SUCCESS);
4066: }
4068: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
4069: /*@
4070: VecLockPop - Pops a read-only lock from a vector
4072: Level: deprecated
4074: .seealso: [](ch_vectors), `Vec`, `VecLockReadPop()`
4075: @*/
4076: PetscErrorCode VecLockPop(Vec x)
4077: {
4078: PetscFunctionBegin;
4079: PetscCall(VecLockReadPop(x));
4080: PetscFunctionReturn(PETSC_SUCCESS);
4081: }
4083: #endif