Actual source code: rvector.c
1: /*
2: Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
3: These are the vector functions the user calls.
4: */
5: #include <petsc/private/vecimpl.h>
7: PetscInt VecGetSubVectorSavedStateId = -1;
9: #if PetscDefined(USE_DEBUG)
10: // this is a no-op '0' macro in optimized builds
11: PetscErrorCode VecValidValues_Internal(Vec vec, PetscInt argnum, PetscBool begin)
12: {
13: PetscFunctionBegin;
14: if (vec->petscnative || vec->ops->getarray) {
15: PetscInt n;
16: const PetscScalar *x;
17: PetscOffloadMask mask;
19: PetscCall(VecGetOffloadMask(vec, &mask));
20: if (!PetscOffloadHost(mask)) PetscFunctionReturn(PETSC_SUCCESS);
21: PetscCall(VecGetLocalSize(vec, &n));
22: PetscCall(VecGetArrayRead(vec, &x));
23: for (PetscInt i = 0; i < n; i++) {
24: if (begin) {
25: PetscCheck(!PetscIsInfOrNanScalar(x[i]), PETSC_COMM_SELF, PETSC_ERR_FP, "Vec entry at local location %" PetscInt_FMT " is not-a-number or infinite at beginning of function: Parameter number %" PetscInt_FMT, i, argnum);
26: } else {
27: PetscCheck(!PetscIsInfOrNanScalar(x[i]), PETSC_COMM_SELF, PETSC_ERR_FP, "Vec entry at local location %" PetscInt_FMT " is not-a-number or infinite at end of function: Parameter number %" PetscInt_FMT, i, argnum);
28: }
29: }
30: PetscCall(VecRestoreArrayRead(vec, &x));
31: }
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
34: #endif
36: /*@
37: VecMaxPointwiseDivide - Computes the maximum of the componentwise division `max = max_i abs(x[i]/y[i])`.
39: Logically Collective
41: Input Parameters:
42: + x - the numerators
43: - y - the denominators
45: Output Parameter:
46: . max - the result
48: Level: advanced
50: Notes:
51: `x` and `y` may be the same vector
53: if a particular `y[i]` is zero, it is treated as 1 in the above formula
55: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`
56: @*/
57: PetscErrorCode VecMaxPointwiseDivide(Vec x, Vec y, PetscReal *max)
58: {
59: PetscFunctionBegin;
62: PetscAssertPointer(max, 3);
65: PetscCheckSameTypeAndComm(x, 1, y, 2);
66: VecCheckSameSize(x, 1, y, 2);
67: VecCheckAssembled(x);
68: VecCheckAssembled(y);
69: PetscCall(VecLockReadPush(x));
70: PetscCall(VecLockReadPush(y));
71: PetscUseTypeMethod(x, maxpointwisedivide, y, max);
72: PetscCall(VecLockReadPop(x));
73: PetscCall(VecLockReadPop(y));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: /*@
78: VecDot - Computes the vector dot product.
80: Collective
82: Input Parameters:
83: + x - first vector
84: - y - second vector
86: Output Parameter:
87: . val - the dot product
89: Level: intermediate
91: Notes for Users of Complex Numbers:
92: For complex vectors, `VecDot()` computes
93: $ val = (x,y) = y^H x,
94: where y^H denotes the conjugate transpose of y. Note that this corresponds to the usual "mathematicians" complex
95: inner product where the SECOND argument gets the complex conjugate. Since the `BLASdot()` complex conjugates the first
96: first argument we call the `BLASdot()` with the arguments reversed.
98: Use `VecTDot()` for the indefinite form
99: $ val = (x,y) = y^T x,
100: where y^T denotes the transpose of y.
102: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDotRealPart()`
103: @*/
104: PetscErrorCode VecDot(Vec x, Vec y, PetscScalar *val)
105: {
106: PetscFunctionBegin;
109: PetscAssertPointer(val, 3);
112: PetscCheckSameTypeAndComm(x, 1, y, 2);
113: VecCheckSameSize(x, 1, y, 2);
114: VecCheckAssembled(x);
115: VecCheckAssembled(y);
117: PetscCall(VecLockReadPush(x));
118: PetscCall(VecLockReadPush(y));
119: PetscCall(PetscLogEventBegin(VEC_Dot, x, y, 0, 0));
120: PetscUseTypeMethod(x, dot, y, val);
121: PetscCall(PetscLogEventEnd(VEC_Dot, x, y, 0, 0));
122: PetscCall(VecLockReadPop(x));
123: PetscCall(VecLockReadPop(y));
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: /*@
128: VecDotRealPart - Computes the real part of the vector dot product.
130: Collective
132: Input Parameters:
133: + x - first vector
134: - y - second vector
136: Output Parameter:
137: . val - the real part of the dot product;
139: Level: intermediate
141: Notes for Users of Complex Numbers:
142: See `VecDot()` for more details on the definition of the dot product for complex numbers
144: For real numbers this returns the same value as `VecDot()`
146: 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
147: the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)
149: Developer Notes:
150: This is not currently optimized to compute only the real part of the dot product.
152: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDot()`, `VecDotNorm2()`
153: @*/
154: PetscErrorCode VecDotRealPart(Vec x, Vec y, PetscReal *val)
155: {
156: PetscScalar fdot;
158: PetscFunctionBegin;
159: PetscCall(VecDot(x, y, &fdot));
160: *val = PetscRealPart(fdot);
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: /*@
165: VecNorm - Computes the vector norm.
167: Collective
169: Input Parameters:
170: + x - the vector
171: - type - the type of the norm requested
173: Output Parameter:
174: . val - the norm
176: Level: intermediate
178: Notes:
179: See `NormType` for descriptions of each norm.
181: For complex numbers `NORM_1` will return the traditional 1 norm of the 2 norm of the complex
182: numbers; that is the 1 norm of the absolute values of the complex entries. In PETSc 3.6 and
183: earlier releases it returned the 1 norm of the 1 norm of the complex entries (what is
184: returned by the BLAS routine `asum()`). Both are valid norms but most people expect the former.
186: This routine stashes the computed norm value, repeated calls before the vector entries are
187: changed are then rapid since the precomputed value is immediately available. Certain vector
188: operations such as `VecSet()` store the norms so the value is immediately available and does
189: not need to be explicitly computed. `VecScale()` updates any stashed norm values, thus calls
190: after `VecScale()` do not need to explicitly recompute the norm.
192: .seealso: [](ch_vectors), `Vec`, `NormType`, `VecDot()`, `VecTDot()`, `VecDotBegin()`, `VecDotEnd()`, `VecNormAvailable()`,
193: `VecNormBegin()`, `VecNormEnd()`, `NormType()`
194: @*/
195: PetscErrorCode VecNorm(Vec x, NormType type, PetscReal *val)
196: {
197: PetscBool flg = PETSC_TRUE;
199: PetscFunctionBegin;
202: VecCheckAssembled(x);
204: PetscAssertPointer(val, 3);
206: PetscCall(VecNormAvailable(x, type, &flg, val));
207: // check that all MPI processes call this routine together and have same availability
208: if (PetscDefined(USE_DEBUG)) {
209: PetscMPIInt b0 = (PetscMPIInt)flg, b1[2], b2[2];
210: b1[0] = -b0;
211: b1[1] = b0;
212: PetscCall(MPIU_Allreduce(b1, b2, 2, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
213: PetscCheck(-b2[0] == b2[1], PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_WRONGSTATE, "Some MPI processes have cached %s norm, others do not. This may happen when some MPI processes call VecGetArray() and some others do not.", NormTypes[type]);
214: if (flg) {
215: PetscReal b1[2], b2[2];
216: b1[0] = -(*val);
217: b1[1] = *val;
218: PetscCall(MPIU_Allreduce(b1, b2, 2, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)x)));
219: PetscCheck(-b2[0] == b2[1], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Difference in cached %s norms: local %g", NormTypes[type], (double)*val);
220: }
221: }
222: if (flg) PetscFunctionReturn(PETSC_SUCCESS);
224: PetscCall(VecLockReadPush(x));
225: PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
226: PetscUseTypeMethod(x, norm, type, val);
227: PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
228: PetscCall(VecLockReadPop(x));
230: if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: /*@
235: VecNormAvailable - Returns the vector norm if it is already known. That is, it has been previously computed and cached in the vector
237: Not Collective
239: Input Parameters:
240: + x - the vector
241: - 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
242: `NORM_1_AND_2`, which computes both norms and stores them
243: in a two element array.
245: Output Parameters:
246: + available - `PETSC_TRUE` if the val returned is valid
247: - val - the norm
249: Level: intermediate
251: Developer Notes:
252: `PETSC_HAVE_SLOW_BLAS_NORM2` will cause a C (loop unrolled) version of the norm to be used, rather
253: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
254: (as opposed to vendor provided) because the FORTRAN BLAS `NRM2()` routine is very slow.
256: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
257: `VecNormBegin()`, `VecNormEnd()`
258: @*/
259: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
260: {
261: PetscFunctionBegin;
264: PetscAssertPointer(available, 3);
265: PetscAssertPointer(val, 4);
267: if (type == NORM_1_AND_2) {
268: *available = PETSC_FALSE;
269: } else {
270: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
271: }
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: /*@
276: VecNormalize - Normalizes a vector by its 2-norm.
278: Collective
280: Input Parameter:
281: . x - the vector
283: Output Parameter:
284: . val - the vector norm before normalization. May be `NULL` if the value is not needed.
286: Level: intermediate
288: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
289: @*/
290: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
291: {
292: PetscReal norm;
294: PetscFunctionBegin;
297: PetscCall(VecSetErrorIfLocked(x, 1));
298: if (val) PetscAssertPointer(val, 2);
299: PetscCall(PetscLogEventBegin(VEC_Normalize, x, 0, 0, 0));
300: PetscCall(VecNorm(x, NORM_2, &norm));
301: if (norm == 0.0) PetscCall(PetscInfo(x, "Vector of zero norm can not be normalized; Returning only the zero norm\n"));
302: else if (PetscIsInfOrNanReal(norm)) PetscCall(PetscInfo(x, "Vector with Inf or Nan norm can not be normalized; Returning only the norm\n"));
303: else {
304: PetscScalar s = 1.0 / norm;
305: PetscCall(VecScale(x, s));
306: }
307: PetscCall(PetscLogEventEnd(VEC_Normalize, x, 0, 0, 0));
308: if (val) *val = norm;
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: /*@C
313: VecMax - Determines the vector component with maximum real part and its location.
315: Collective
317: Input Parameter:
318: . x - the vector
320: Output Parameters:
321: + p - the index of `val` (pass `NULL` if you don't want this) in the vector
322: - val - the maximum component
324: Level: intermediate
326: Notes:
327: Returns the value `PETSC_MIN_REAL` and negative `p` if the vector is of length 0.
329: Returns the smallest index with the maximum value
331: Developer Note:
332: The Nag Fortran compiler does not like the symbol name VecMax
334: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
335: @*/
336: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
337: {
338: PetscFunctionBegin;
341: VecCheckAssembled(x);
342: if (p) PetscAssertPointer(p, 2);
343: PetscAssertPointer(val, 3);
344: PetscCall(VecLockReadPush(x));
345: PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
346: PetscUseTypeMethod(x, max, p, val);
347: PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
348: PetscCall(VecLockReadPop(x));
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: /*@C
353: VecMin - Determines the vector component with minimum real part and its location.
355: Collective
357: Input Parameter:
358: . x - the vector
360: Output Parameters:
361: + p - the index of `val` (pass `NULL` if you don't want this location) in the vector
362: - val - the minimum component
364: Level: intermediate
366: Notes:
367: Returns the value `PETSC_MAX_REAL` and negative `p` if the vector is of length 0.
369: This returns the smallest index with the minimum value
371: Developer Note:
372: The Nag Fortran compiler does not like the symbol name VecMin
374: .seealso: [](ch_vectors), `Vec`, `VecMax()`
375: @*/
376: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
377: {
378: PetscFunctionBegin;
381: VecCheckAssembled(x);
382: if (p) PetscAssertPointer(p, 2);
383: PetscAssertPointer(val, 3);
384: PetscCall(VecLockReadPush(x));
385: PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
386: PetscUseTypeMethod(x, min, p, val);
387: PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
388: PetscCall(VecLockReadPop(x));
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: /*@
393: VecTDot - Computes an indefinite vector dot product. That is, this
394: routine does NOT use the complex conjugate.
396: Collective
398: Input Parameters:
399: + x - first vector
400: - y - second vector
402: Output Parameter:
403: . val - the dot product
405: Level: intermediate
407: Notes for Users of Complex Numbers:
408: For complex vectors, `VecTDot()` computes the indefinite form
409: $ val = (x,y) = y^T x,
410: where y^T denotes the transpose of y.
412: Use `VecDot()` for the inner product
413: $ val = (x,y) = y^H x,
414: where y^H denotes the conjugate transpose of y.
416: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
417: @*/
418: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
419: {
420: PetscFunctionBegin;
423: PetscAssertPointer(val, 3);
426: PetscCheckSameTypeAndComm(x, 1, y, 2);
427: VecCheckSameSize(x, 1, y, 2);
428: VecCheckAssembled(x);
429: VecCheckAssembled(y);
431: PetscCall(VecLockReadPush(x));
432: PetscCall(VecLockReadPush(y));
433: PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
434: PetscUseTypeMethod(x, tdot, y, val);
435: PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
436: PetscCall(VecLockReadPop(x));
437: PetscCall(VecLockReadPop(y));
438: PetscFunctionReturn(PETSC_SUCCESS);
439: }
441: PetscErrorCode VecScaleAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
442: {
443: PetscReal norms[4];
444: PetscBool flgs[4];
445: PetscScalar one = 1.0;
447: PetscFunctionBegin;
450: VecCheckAssembled(x);
451: PetscCall(VecSetErrorIfLocked(x, 1));
453: if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);
455: /* get current stashed norms */
456: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
458: PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
459: VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
460: PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));
462: PetscCall(PetscObjectStateIncrease((PetscObject)x));
463: /* put the scaled stashed norms back into the Vec */
464: for (PetscInt i = 0; i < 4; i++) {
465: PetscReal ar = PetscAbsScalar(alpha);
466: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
467: }
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
471: /*@
472: VecScale - Scales a vector.
474: Logically Collective
476: Input Parameters:
477: + x - the vector
478: - alpha - the scalar
480: Level: intermediate
482: Note:
483: For a vector with n components, `VecScale()` computes x[i] = alpha * x[i], for i=1,...,n.
485: .seealso: [](ch_vectors), `Vec`, `VecSet()`
486: @*/
487: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
488: {
489: PetscFunctionBegin;
490: PetscCall(VecScaleAsync_Private(x, alpha, NULL));
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
495: {
496: PetscFunctionBegin;
499: VecCheckAssembled(x);
501: PetscCall(VecSetErrorIfLocked(x, 1));
503: if (alpha == 0) {
504: PetscReal norm;
505: PetscBool set;
507: PetscCall(VecNormAvailable(x, NORM_2, &set, &norm));
508: if (set == PETSC_TRUE && norm == 0) PetscFunctionReturn(PETSC_SUCCESS);
509: }
510: PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
511: VecMethodDispatch(x, dctx, VecAsyncFnName(Set), set, (Vec, PetscScalar, PetscDeviceContext), alpha);
512: PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
513: PetscCall(PetscObjectStateIncrease((PetscObject)x));
515: /* norms can be simply set (if |alpha|*N not too large) */
516: {
517: PetscReal val = PetscAbsScalar(alpha);
518: const PetscInt N = x->map->N;
520: if (N == 0) {
521: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0));
522: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
523: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
524: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
525: } else if (val > PETSC_MAX_REAL / N) {
526: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
527: } else {
528: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
529: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
530: val *= PetscSqrtReal((PetscReal)N);
531: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
532: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
533: }
534: }
535: PetscFunctionReturn(PETSC_SUCCESS);
536: }
538: /*@
539: VecSet - Sets all components of a vector to a single scalar value.
541: Logically Collective
543: Input Parameters:
544: + x - the vector
545: - alpha - the scalar
547: Level: beginner
549: Notes:
550: For a vector of dimension n, `VecSet()` sets x[i] = alpha, for i=1,...,n,
551: so that all vector entries then equal the identical
552: scalar value, `alpha`. Use the more general routine
553: `VecSetValues()` to set different vector entries.
555: You CANNOT call this after you have called `VecSetValues()` but before you call
556: `VecAssemblyBegin()`
558: If `alpha` is zero and the norm of the vector is known to be zero then this skips the unneeded zeroing process
560: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
561: @*/
562: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
563: {
564: PetscFunctionBegin;
565: PetscCall(VecSetAsync_Private(x, alpha, NULL));
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: PetscErrorCode VecAXPYAsync_Private(Vec y, PetscScalar alpha, Vec x, PetscDeviceContext dctx)
570: {
571: PetscFunctionBegin;
576: PetscCheckSameTypeAndComm(x, 3, y, 1);
577: VecCheckSameSize(x, 3, y, 1);
578: VecCheckAssembled(x);
579: VecCheckAssembled(y);
581: if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
582: PetscCall(VecSetErrorIfLocked(y, 1));
583: if (x == y) {
584: PetscCall(VecScale(y, alpha + 1.0));
585: PetscFunctionReturn(PETSC_SUCCESS);
586: }
587: PetscCall(VecLockReadPush(x));
588: PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
589: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPY), axpy, (Vec, PetscScalar, Vec, PetscDeviceContext), alpha, x);
590: PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
591: PetscCall(VecLockReadPop(x));
592: PetscCall(PetscObjectStateIncrease((PetscObject)y));
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
595: /*@
596: VecAXPY - Computes `y = alpha x + y`.
598: Logically Collective
600: Input Parameters:
601: + alpha - the scalar
602: . x - vector scale by `alpha`
603: - y - vector accumulated into
605: Output Parameter:
606: . y - output vector
608: Level: intermediate
610: Notes:
611: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
612: .vb
613: VecAXPY(y,alpha,x) y = alpha x + y
614: VecAYPX(y,beta,x) y = x + beta y
615: VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
616: VecWAXPY(w,alpha,x,y) w = alpha x + y
617: VecAXPBYPCZ(z,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
618: VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
619: .ve
621: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
622: @*/
623: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
624: {
625: PetscFunctionBegin;
626: PetscCall(VecAXPYAsync_Private(y, alpha, x, NULL));
627: PetscFunctionReturn(PETSC_SUCCESS);
628: }
630: PetscErrorCode VecAYPXAsync_Private(Vec y, PetscScalar beta, Vec x, PetscDeviceContext dctx)
631: {
632: PetscFunctionBegin;
637: PetscCheckSameTypeAndComm(x, 3, y, 1);
638: VecCheckSameSize(x, 1, y, 3);
639: VecCheckAssembled(x);
640: VecCheckAssembled(y);
642: PetscCall(VecSetErrorIfLocked(y, 1));
643: if (x == y) {
644: PetscCall(VecScale(y, beta + 1.0));
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
647: PetscCall(VecLockReadPush(x));
648: if (beta == (PetscScalar)0.0) {
649: PetscCall(VecCopy(x, y));
650: } else {
651: PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
652: VecMethodDispatch(y, dctx, VecAsyncFnName(AYPX), aypx, (Vec, PetscScalar, Vec, PetscDeviceContext), beta, x);
653: PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
654: PetscCall(PetscObjectStateIncrease((PetscObject)y));
655: }
656: PetscCall(VecLockReadPop(x));
657: PetscFunctionReturn(PETSC_SUCCESS);
658: }
660: /*@
661: VecAYPX - Computes `y = x + beta y`.
663: Logically Collective
665: Input Parameters:
666: + beta - the scalar
667: . x - the unscaled vector
668: - y - the vector to be scaled
670: Output Parameter:
671: . y - output vector
673: Level: intermediate
675: Developer Notes:
676: The implementation is optimized for `beta` of -1.0, 0.0, and 1.0
678: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
679: @*/
680: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
681: {
682: PetscFunctionBegin;
683: PetscCall(VecAYPXAsync_Private(y, beta, x, NULL));
684: PetscFunctionReturn(PETSC_SUCCESS);
685: }
687: PetscErrorCode VecAXPBYAsync_Private(Vec y, PetscScalar alpha, PetscScalar beta, Vec x, PetscDeviceContext dctx)
688: {
689: PetscFunctionBegin;
694: PetscCheckSameTypeAndComm(x, 4, y, 1);
695: VecCheckSameSize(y, 1, x, 4);
696: VecCheckAssembled(x);
697: VecCheckAssembled(y);
700: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
701: if (x == y) {
702: PetscCall(VecScale(y, alpha + beta));
703: PetscFunctionReturn(PETSC_SUCCESS);
704: }
706: PetscCall(VecSetErrorIfLocked(y, 1));
707: PetscCall(VecLockReadPush(x));
708: PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
709: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPBY), axpby, (Vec, PetscScalar, PetscScalar, Vec, PetscDeviceContext), alpha, beta, x);
710: PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
711: PetscCall(PetscObjectStateIncrease((PetscObject)y));
712: PetscCall(VecLockReadPop(x));
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }
716: /*@
717: VecAXPBY - Computes `y = alpha x + beta y`.
719: Logically Collective
721: Input Parameters:
722: + alpha - first scalar
723: . beta - second scalar
724: . x - the first scaled vector
725: - y - the second scaled vector
727: Output Parameter:
728: . y - output vector
730: Level: intermediate
732: Developer Notes:
733: The implementation is optimized for `alpha` and/or `beta` values of 0.0 and 1.0
735: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
736: @*/
737: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
738: {
739: PetscFunctionBegin;
740: PetscCall(VecAXPBYAsync_Private(y, alpha, beta, x, NULL));
741: PetscFunctionReturn(PETSC_SUCCESS);
742: }
744: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
745: {
746: PetscFunctionBegin;
753: PetscCheckSameTypeAndComm(x, 5, y, 6);
754: PetscCheckSameTypeAndComm(x, 5, z, 1);
755: VecCheckSameSize(x, 5, y, 6);
756: VecCheckSameSize(x, 5, z, 1);
757: PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
758: PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
759: VecCheckAssembled(x);
760: VecCheckAssembled(y);
761: VecCheckAssembled(z);
765: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
767: PetscCall(VecSetErrorIfLocked(z, 1));
768: PetscCall(VecLockReadPush(x));
769: PetscCall(VecLockReadPush(y));
770: PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
771: VecMethodDispatch(z, dctx, VecAsyncFnName(AXPBYPCZ), axpbypcz, (Vec, PetscScalar, PetscScalar, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, beta, gamma, x, y);
772: PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
773: PetscCall(PetscObjectStateIncrease((PetscObject)z));
774: PetscCall(VecLockReadPop(x));
775: PetscCall(VecLockReadPop(y));
776: PetscFunctionReturn(PETSC_SUCCESS);
777: }
778: /*@
779: VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`
781: Logically Collective
783: Input Parameters:
784: + alpha - first scalar
785: . beta - second scalar
786: . gamma - third scalar
787: . x - first vector
788: . y - second vector
789: - z - third vector
791: Output Parameter:
792: . z - output vector
794: Level: intermediate
796: Note:
797: `x`, `y` and `z` must be different vectors
799: Developer Notes:
800: The implementation is optimized for `alpha` of 1.0 and `gamma` of 1.0 or 0.0
802: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
803: @*/
804: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
805: {
806: PetscFunctionBegin;
807: PetscCall(VecAXPBYPCZAsync_Private(z, alpha, beta, gamma, x, y, NULL));
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
812: {
813: PetscFunctionBegin;
820: PetscCheckSameTypeAndComm(x, 3, y, 4);
821: PetscCheckSameTypeAndComm(y, 4, w, 1);
822: VecCheckSameSize(x, 3, y, 4);
823: VecCheckSameSize(x, 3, w, 1);
824: PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
825: PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
826: VecCheckAssembled(x);
827: VecCheckAssembled(y);
829: PetscCall(VecSetErrorIfLocked(w, 1));
831: PetscCall(VecLockReadPush(x));
832: PetscCall(VecLockReadPush(y));
833: if (alpha == (PetscScalar)0.0) {
834: PetscCall(VecCopyAsync_Private(y, w, dctx));
835: } else {
836: PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
837: VecMethodDispatch(w, dctx, VecAsyncFnName(WAXPY), waxpy, (Vec, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, x, y);
838: PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
839: PetscCall(PetscObjectStateIncrease((PetscObject)w));
840: }
841: PetscCall(VecLockReadPop(x));
842: PetscCall(VecLockReadPop(y));
843: PetscFunctionReturn(PETSC_SUCCESS);
844: }
846: /*@
847: VecWAXPY - Computes `w = alpha x + y`.
849: Logically Collective
851: Input Parameters:
852: + alpha - the scalar
853: . x - first vector, multiplied by `alpha`
854: - y - second vector
856: Output Parameter:
857: . w - the result
859: Level: intermediate
861: Note:
862: `w` cannot be either `x` or `y`, but `x` and `y` can be the same
864: Developer Notes:
865: The implementation is optimized for alpha of -1.0, 0.0, and 1.0
867: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
868: @*/
869: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
870: {
871: PetscFunctionBegin;
872: PetscCall(VecWAXPYAsync_Private(w, alpha, x, y, NULL));
873: PetscFunctionReturn(PETSC_SUCCESS);
874: }
876: /*@C
877: VecSetValues - Inserts or adds values into certain locations of a vector.
879: Not Collective
881: Input Parameters:
882: + x - vector to insert in
883: . ni - number of elements to add
884: . ix - indices where to add
885: . y - array of values
886: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries
888: Level: beginner
890: Notes:
891: .vb
892: `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
893: .ve
895: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
896: options cannot be mixed without intervening calls to the assembly
897: routines.
899: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
900: MUST be called after all calls to `VecSetValues()` have been completed.
902: VecSetValues() uses 0-based indices in Fortran as well as in C.
904: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
905: negative indices may be passed in ix. These rows are
906: simply ignored. This allows easily inserting element load matrices
907: with homogeneous Dirichlet boundary conditions that you don't want represented
908: in the vector.
910: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
911: `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
912: @*/
913: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
914: {
915: PetscFunctionBeginHot;
917: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
918: PetscAssertPointer(ix, 3);
919: PetscAssertPointer(y, 4);
922: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
923: PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
924: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
925: PetscCall(PetscObjectStateIncrease((PetscObject)x));
926: PetscFunctionReturn(PETSC_SUCCESS);
927: }
929: /*@C
930: VecGetValues - Gets values from certain locations of a vector. Currently
931: can only get values on the same processor on which they are owned
933: Not Collective
935: Input Parameters:
936: + x - vector to get values from
937: . ni - number of elements to get
938: - ix - indices where to get them from (in global 1d numbering)
940: Output Parameter:
941: . y - array of values
943: Level: beginner
945: Notes:
946: The user provides the allocated array y; it is NOT allocated in this routine
948: `VecGetValues()` gets y[i] = x[ix[i]], for i=0,...,ni-1.
950: `VecAssemblyBegin()` and `VecAssemblyEnd()` MUST be called before calling this if `VecSetValues()` or related routine has been called
952: VecGetValues() uses 0-based indices in Fortran as well as in C.
954: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
955: negative indices may be passed in ix. These rows are
956: simply ignored.
958: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
959: @*/
960: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
961: {
962: PetscFunctionBegin;
964: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
965: PetscAssertPointer(ix, 3);
966: PetscAssertPointer(y, 4);
968: VecCheckAssembled(x);
969: PetscUseTypeMethod(x, getvalues, ni, ix, y);
970: PetscFunctionReturn(PETSC_SUCCESS);
971: }
973: /*@C
974: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
976: Not Collective
978: Input Parameters:
979: + x - vector to insert in
980: . ni - number of blocks to add
981: . ix - indices where to add in block count, rather than element count
982: . y - array of values
983: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries
985: Level: intermediate
987: Notes:
988: `VecSetValuesBlocked()` sets x[bs*ix[i]+j] = y[bs*i+j],
989: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
991: Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
992: options cannot be mixed without intervening calls to the assembly
993: routines.
995: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
996: MUST be called after all calls to `VecSetValuesBlocked()` have been completed.
998: `VecSetValuesBlocked()` uses 0-based indices in Fortran as well as in C.
1000: Negative indices may be passed in ix, these rows are
1001: simply ignored. This allows easily inserting element load matrices
1002: with homogeneous Dirichlet boundary conditions that you don't want represented
1003: in the vector.
1005: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
1006: `VecSetValues()`
1007: @*/
1008: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1009: {
1010: PetscFunctionBeginHot;
1012: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1013: PetscAssertPointer(ix, 3);
1014: PetscAssertPointer(y, 4);
1017: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1018: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1019: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1020: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1021: PetscFunctionReturn(PETSC_SUCCESS);
1022: }
1024: /*@C
1025: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1026: using a local ordering of the nodes.
1028: Not Collective
1030: Input Parameters:
1031: + x - vector to insert in
1032: . ni - number of elements to add
1033: . ix - indices where to add
1034: . y - array of values
1035: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1037: Level: intermediate
1039: Notes:
1040: `VecSetValuesLocal()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
1042: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
1043: options cannot be mixed without intervening calls to the assembly
1044: routines.
1046: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1047: MUST be called after all calls to `VecSetValuesLocal()` have been completed.
1049: `VecSetValuesLocal()` uses 0-based indices in Fortran as well as in C.
1051: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1052: `VecSetValuesBlockedLocal()`
1053: @*/
1054: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1055: {
1056: PetscInt lixp[128], *lix = lixp;
1058: PetscFunctionBeginHot;
1060: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1061: PetscAssertPointer(ix, 3);
1062: PetscAssertPointer(y, 4);
1065: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1066: if (!x->ops->setvalueslocal) {
1067: if (x->map->mapping) {
1068: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1069: PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1070: PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1071: if (ni > 128) PetscCall(PetscFree(lix));
1072: } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1073: } else PetscUseTypeMethod(x, setvalueslocal, ni, ix, y, iora);
1074: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1075: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1076: PetscFunctionReturn(PETSC_SUCCESS);
1077: }
1079: /*@
1080: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1081: using a local ordering of the nodes.
1083: Not Collective
1085: Input Parameters:
1086: + x - vector to insert in
1087: . ni - number of blocks to add
1088: . ix - indices where to add in block count, not element count
1089: . y - array of values
1090: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1092: Level: intermediate
1094: Notes:
1095: `VecSetValuesBlockedLocal()` sets x[bs*ix[i]+j] = y[bs*i+j],
1096: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with `VecSetBlockSize()`.
1098: Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1099: options cannot be mixed without intervening calls to the assembly
1100: routines.
1102: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1103: MUST be called after all calls to `VecSetValuesBlockedLocal()` have been completed.
1105: `VecSetValuesBlockedLocal()` uses 0-based indices in Fortran as well as in C.
1107: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1108: `VecSetLocalToGlobalMapping()`
1109: @*/
1110: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1111: {
1112: PetscInt lixp[128], *lix = lixp;
1114: PetscFunctionBeginHot;
1116: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1117: PetscAssertPointer(ix, 3);
1118: PetscAssertPointer(y, 4);
1120: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1121: if (x->map->mapping) {
1122: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1123: PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1124: PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1125: if (ni > 128) PetscCall(PetscFree(lix));
1126: } else {
1127: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1128: }
1129: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1130: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1131: PetscFunctionReturn(PETSC_SUCCESS);
1132: }
1134: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1135: {
1136: PetscFunctionBegin;
1139: VecCheckAssembled(x);
1141: if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1142: PetscAssertPointer(y, 3);
1143: for (PetscInt i = 0; i < nv; ++i) {
1146: PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1147: VecCheckSameSize(x, 1, y[i], 3);
1148: VecCheckAssembled(y[i]);
1149: PetscCall(VecLockReadPush(y[i]));
1150: }
1151: PetscAssertPointer(result, 4);
1154: PetscCall(VecLockReadPush(x));
1155: PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1156: PetscCall((*mxdot)(x, nv, y, result));
1157: PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1158: PetscCall(VecLockReadPop(x));
1159: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1160: PetscFunctionReturn(PETSC_SUCCESS);
1161: }
1163: /*@
1164: VecMTDot - Computes indefinite vector multiple dot products.
1165: That is, it does NOT use the complex conjugate.
1167: Collective
1169: Input Parameters:
1170: + x - one vector
1171: . nv - number of vectors
1172: - y - array of vectors. Note that vectors are pointers
1174: Output Parameter:
1175: . val - array of the dot products
1177: Level: intermediate
1179: Notes for Users of Complex Numbers:
1180: For complex vectors, `VecMTDot()` computes the indefinite form
1181: $ val = (x,y) = y^T x,
1182: where y^T denotes the transpose of y.
1184: Use `VecMDot()` for the inner product
1185: $ val = (x,y) = y^H x,
1186: where y^H denotes the conjugate transpose of y.
1188: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1189: @*/
1190: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1191: {
1192: PetscFunctionBegin;
1194: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1195: PetscFunctionReturn(PETSC_SUCCESS);
1196: }
1198: /*@
1199: VecMDot - Computes multiple vector dot products.
1201: Collective
1203: Input Parameters:
1204: + x - one vector
1205: . nv - number of vectors
1206: - y - array of vectors.
1208: Output Parameter:
1209: . val - array of the dot products (does not allocate the array)
1211: Level: intermediate
1213: Notes for Users of Complex Numbers:
1214: For complex vectors, `VecMDot()` computes
1215: $ val = (x,y) = y^H x,
1216: where y^H denotes the conjugate transpose of y.
1218: Use `VecMTDot()` for the indefinite form
1219: $ val = (x,y) = y^T x,
1220: where y^T denotes the transpose of y.
1222: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`
1223: @*/
1224: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1225: {
1226: PetscFunctionBegin;
1228: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1229: PetscFunctionReturn(PETSC_SUCCESS);
1230: }
1232: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
1233: {
1234: PetscFunctionBegin;
1236: VecCheckAssembled(y);
1238: PetscCall(VecSetErrorIfLocked(y, 1));
1239: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1240: if (nv) {
1241: PetscInt zeros = 0;
1243: PetscAssertPointer(alpha, 3);
1244: PetscAssertPointer(x, 4);
1245: for (PetscInt i = 0; i < nv; ++i) {
1249: PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1250: VecCheckSameSize(y, 1, x[i], 4);
1251: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1252: VecCheckAssembled(x[i]);
1253: PetscCall(VecLockReadPush(x[i]));
1254: zeros += alpha[i] == (PetscScalar)0.0;
1255: }
1257: if (zeros < nv) {
1258: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1259: VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1260: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1261: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1262: }
1264: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1265: }
1266: PetscFunctionReturn(PETSC_SUCCESS);
1267: }
1269: /*@
1270: VecMAXPY - Computes `y = y + sum alpha[i] x[i]`
1272: Logically Collective
1274: Input Parameters:
1275: + nv - number of scalars and x-vectors
1276: . alpha - array of scalars
1277: . y - one vector
1278: - x - array of vectors
1280: Level: intermediate
1282: Note:
1283: `y` cannot be any of the `x` vectors
1285: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`,`VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1286: @*/
1287: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1288: {
1289: PetscFunctionBegin;
1290: PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1291: PetscFunctionReturn(PETSC_SUCCESS);
1292: }
1294: /*@
1295: VecMAXPBY - Computes `y = beta y + sum alpha[i] x[i]`
1297: Logically Collective
1299: Input Parameters:
1300: + nv - number of scalars and x-vectors
1301: . alpha - array of scalars
1302: . beta - scalar
1303: . y - one vector
1304: - x - array of vectors
1306: Level: intermediate
1308: Note:
1309: `y` cannot be any of the `x` vectors.
1311: Developer Notes:
1312: This is a convenience routine, but implementations might be able to optimize it, for example, when `beta` is zero.
1314: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1315: @*/
1316: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
1317: {
1318: PetscFunctionBegin;
1320: VecCheckAssembled(y);
1322: PetscCall(VecSetErrorIfLocked(y, 1));
1323: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1326: if (y->ops->maxpby) {
1327: PetscInt zeros = 0;
1329: if (nv) {
1330: PetscAssertPointer(alpha, 3);
1331: PetscAssertPointer(x, 5);
1332: }
1334: for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
1338: PetscCheckSameTypeAndComm(y, 1, x[i], 5);
1339: VecCheckSameSize(y, 1, x[i], 5);
1340: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1341: VecCheckAssembled(x[i]);
1342: PetscCall(VecLockReadPush(x[i]));
1343: zeros += alpha[i] == (PetscScalar)0.0;
1344: }
1346: if (zeros < nv) { // has nonzero alpha
1347: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1348: PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1349: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1350: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1351: } else {
1352: PetscCall(VecScale(y, beta));
1353: }
1355: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1356: } else { // no maxpby
1357: if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1358: else PetscCall(VecScale(y, beta));
1359: PetscCall(VecMAXPY(y, nv, alpha, x));
1360: }
1361: PetscFunctionReturn(PETSC_SUCCESS);
1362: }
1364: /*@
1365: VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
1366: in the order they appear in the array. The concatenated vector resides on the same
1367: communicator and is the same type as the source vectors.
1369: Collective
1371: Input Parameters:
1372: + nx - number of vectors to be concatenated
1373: - X - array containing the vectors to be concatenated in the order of concatenation
1375: Output Parameters:
1376: + Y - concatenated vector
1377: - x_is - array of index sets corresponding to the concatenated components of `Y` (pass `NULL` if not needed)
1379: Level: advanced
1381: Notes:
1382: Concatenation is similar to the functionality of a `VECNEST` object; they both represent combination of
1383: different vector spaces. However, concatenated vectors do not store any information about their
1384: sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1385: manipulation of data in the concatenated vector that corresponds to the original components at creation.
1387: This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1388: has to operate on combined vector spaces and cannot utilize `VECNEST` objects due to incompatibility with
1389: bound projections.
1391: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1392: @*/
1393: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1394: {
1395: MPI_Comm comm;
1396: VecType vec_type;
1397: Vec Ytmp, Xtmp;
1398: IS *is_tmp;
1399: PetscInt i, shift = 0, Xnl, Xng, Xbegin;
1401: PetscFunctionBegin;
1405: PetscAssertPointer(Y, 3);
1407: if ((*X)->ops->concatenate) {
1408: /* use the dedicated concatenation function if available */
1409: PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1410: } else {
1411: /* loop over vectors and start creating IS */
1412: comm = PetscObjectComm((PetscObject)*X);
1413: PetscCall(VecGetType(*X, &vec_type));
1414: PetscCall(PetscMalloc1(nx, &is_tmp));
1415: for (i = 0; i < nx; i++) {
1416: PetscCall(VecGetSize(X[i], &Xng));
1417: PetscCall(VecGetLocalSize(X[i], &Xnl));
1418: PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1419: PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1420: shift += Xng;
1421: }
1422: /* create the concatenated vector */
1423: PetscCall(VecCreate(comm, &Ytmp));
1424: PetscCall(VecSetType(Ytmp, vec_type));
1425: PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1426: PetscCall(VecSetUp(Ytmp));
1427: /* copy data from X array to Y and return */
1428: for (i = 0; i < nx; i++) {
1429: PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1430: PetscCall(VecCopy(X[i], Xtmp));
1431: PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1432: }
1433: *Y = Ytmp;
1434: if (x_is) {
1435: *x_is = is_tmp;
1436: } else {
1437: for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1438: PetscCall(PetscFree(is_tmp));
1439: }
1440: }
1441: PetscFunctionReturn(PETSC_SUCCESS);
1442: }
1444: /* A helper function for VecGetSubVector to check if we can implement it with no-copy (i.e. the subvector shares
1445: memory with the original vector), and the block size of the subvector.
1447: Input Parameters:
1448: + X - the original vector
1449: - is - the index set of the subvector
1451: Output Parameters:
1452: + contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
1453: . start - start of contiguous block, as an offset from the start of the ownership range of the original vector
1454: - blocksize - the block size of the subvector
1456: */
1457: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1458: {
1459: PetscInt gstart, gend, lstart;
1460: PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1461: PetscInt n, N, ibs, vbs, bs = -1;
1463: PetscFunctionBegin;
1464: PetscCall(ISGetLocalSize(is, &n));
1465: PetscCall(ISGetSize(is, &N));
1466: PetscCall(ISGetBlockSize(is, &ibs));
1467: PetscCall(VecGetBlockSize(X, &vbs));
1468: PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1469: PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1470: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1471: if (ibs > 1) {
1472: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1473: bs = ibs;
1474: } else {
1475: if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1476: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1477: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1478: }
1480: *contig = red[0];
1481: *start = lstart;
1482: *blocksize = bs;
1483: PetscFunctionReturn(PETSC_SUCCESS);
1484: }
1486: /* A helper function for VecGetSubVector, to be used when we have to build a standalone subvector through VecScatter
1488: Input Parameters:
1489: + X - the original vector
1490: . is - the index set of the subvector
1491: - bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()
1493: Output Parameter:
1494: . Z - the subvector, which will compose the VecScatter context on output
1495: */
1496: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1497: {
1498: PetscInt n, N;
1499: VecScatter vscat;
1500: Vec Y;
1502: PetscFunctionBegin;
1503: PetscCall(ISGetLocalSize(is, &n));
1504: PetscCall(ISGetSize(is, &N));
1505: PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1506: PetscCall(VecSetSizes(Y, n, N));
1507: PetscCall(VecSetBlockSize(Y, bs));
1508: PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1509: PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1510: PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1511: PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1512: PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1513: PetscCall(VecScatterDestroy(&vscat));
1514: *Z = Y;
1515: PetscFunctionReturn(PETSC_SUCCESS);
1516: }
1518: /*@
1519: VecGetSubVector - Gets a vector representing part of another vector
1521: Collective
1523: Input Parameters:
1524: + X - vector from which to extract a subvector
1525: - is - index set representing portion of `X` to extract
1527: Output Parameter:
1528: . Y - subvector corresponding to `is`
1530: Level: advanced
1532: Notes:
1533: The subvector `Y` should be returned with `VecRestoreSubVector()`.
1534: `X` and `is` must be defined on the same communicator
1536: Changes to the subvector will be reflected in the `X` vector on the call to `VecRestoreSubVector()`.
1538: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1539: modifying the subvector. Other non-overlapping subvectors can still be obtained from `X` using this function.
1541: 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`.
1543: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1544: @*/
1545: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1546: {
1547: Vec Z;
1549: PetscFunctionBegin;
1552: PetscCheckSameComm(X, 1, is, 2);
1553: PetscAssertPointer(Y, 3);
1554: if (X->ops->getsubvector) {
1555: PetscUseTypeMethod(X, getsubvector, is, &Z);
1556: } else { /* Default implementation currently does no caching */
1557: PetscBool contig;
1558: PetscInt n, N, start, bs;
1560: PetscCall(ISGetLocalSize(is, &n));
1561: PetscCall(ISGetSize(is, &N));
1562: PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1563: if (contig) { /* We can do a no-copy implementation */
1564: const PetscScalar *x;
1565: PetscInt state = 0;
1566: PetscBool isstd, iscuda, iship;
1568: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1569: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1570: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1571: if (iscuda) {
1572: #if defined(PETSC_HAVE_CUDA)
1573: const PetscScalar *x_d;
1574: PetscMPIInt size;
1575: PetscOffloadMask flg;
1577: PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1578: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1579: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1580: if (x) x += start;
1581: if (x_d) x_d += start;
1582: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1583: if (size == 1) {
1584: PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1585: } else {
1586: PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1587: }
1588: Z->offloadmask = flg;
1589: #endif
1590: } else if (iship) {
1591: #if defined(PETSC_HAVE_HIP)
1592: const PetscScalar *x_d;
1593: PetscMPIInt size;
1594: PetscOffloadMask flg;
1596: PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1597: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1598: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1599: if (x) x += start;
1600: if (x_d) x_d += start;
1601: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1602: if (size == 1) {
1603: PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1604: } else {
1605: PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1606: }
1607: Z->offloadmask = flg;
1608: #endif
1609: } else if (isstd) {
1610: PetscMPIInt size;
1612: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1613: PetscCall(VecGetArrayRead(X, &x));
1614: if (x) x += start;
1615: if (size == 1) {
1616: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1617: } else {
1618: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1619: }
1620: PetscCall(VecRestoreArrayRead(X, &x));
1621: } else { /* default implementation: use place array */
1622: PetscCall(VecGetArrayRead(X, &x));
1623: PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1624: PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1625: PetscCall(VecSetSizes(Z, n, N));
1626: PetscCall(VecSetBlockSize(Z, bs));
1627: PetscCall(VecPlaceArray(Z, PetscSafePointerPlusOffset(x, start)));
1628: PetscCall(VecRestoreArrayRead(X, &x));
1629: }
1631: /* this is relevant only in debug mode */
1632: PetscCall(VecLockGet(X, &state));
1633: if (state) PetscCall(VecLockReadPush(Z));
1634: Z->ops->placearray = NULL;
1635: Z->ops->replacearray = NULL;
1636: } else { /* Have to create a scatter and do a copy */
1637: PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1638: }
1639: }
1640: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1641: if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1642: PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1643: *Y = Z;
1644: PetscFunctionReturn(PETSC_SUCCESS);
1645: }
1647: /*@
1648: VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`
1650: Collective
1652: Input Parameters:
1653: + X - vector from which subvector was obtained
1654: . is - index set representing the subset of `X`
1655: - Y - subvector being restored
1657: Level: advanced
1659: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1660: @*/
1661: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1662: {
1663: PETSC_UNUSED PetscObjectState dummystate = 0;
1664: PetscBool unchanged;
1666: PetscFunctionBegin;
1669: PetscCheckSameComm(X, 1, is, 2);
1670: PetscAssertPointer(Y, 3);
1673: if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1674: else {
1675: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1676: if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1677: VecScatter scatter;
1678: PetscInt state;
1680: PetscCall(VecLockGet(X, &state));
1681: PetscCheck(state == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec X is locked for read-only or read/write access");
1683: PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1684: if (scatter) {
1685: PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1686: PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1687: } else {
1688: PetscBool iscuda, iship;
1689: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1690: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1692: if (iscuda) {
1693: #if defined(PETSC_HAVE_CUDA)
1694: PetscOffloadMask ymask = (*Y)->offloadmask;
1696: /* The offloadmask of X dictates where to move memory
1697: If X GPU data is valid, then move Y data on GPU if needed
1698: Otherwise, move back to the CPU */
1699: switch (X->offloadmask) {
1700: case PETSC_OFFLOAD_BOTH:
1701: if (ymask == PETSC_OFFLOAD_CPU) {
1702: PetscCall(VecCUDAResetArray(*Y));
1703: } else if (ymask == PETSC_OFFLOAD_GPU) {
1704: X->offloadmask = PETSC_OFFLOAD_GPU;
1705: }
1706: break;
1707: case PETSC_OFFLOAD_GPU:
1708: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1709: break;
1710: case PETSC_OFFLOAD_CPU:
1711: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1712: break;
1713: case PETSC_OFFLOAD_UNALLOCATED:
1714: case PETSC_OFFLOAD_KOKKOS:
1715: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1716: }
1717: #endif
1718: } else if (iship) {
1719: #if defined(PETSC_HAVE_HIP)
1720: PetscOffloadMask ymask = (*Y)->offloadmask;
1722: /* The offloadmask of X dictates where to move memory
1723: If X GPU data is valid, then move Y data on GPU if needed
1724: Otherwise, move back to the CPU */
1725: switch (X->offloadmask) {
1726: case PETSC_OFFLOAD_BOTH:
1727: if (ymask == PETSC_OFFLOAD_CPU) {
1728: PetscCall(VecHIPResetArray(*Y));
1729: } else if (ymask == PETSC_OFFLOAD_GPU) {
1730: X->offloadmask = PETSC_OFFLOAD_GPU;
1731: }
1732: break;
1733: case PETSC_OFFLOAD_GPU:
1734: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1735: break;
1736: case PETSC_OFFLOAD_CPU:
1737: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1738: break;
1739: case PETSC_OFFLOAD_UNALLOCATED:
1740: case PETSC_OFFLOAD_KOKKOS:
1741: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1742: }
1743: #endif
1744: } else {
1745: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1746: PetscCall(VecResetArray(*Y));
1747: }
1748: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1749: }
1750: }
1751: }
1752: PetscCall(VecDestroy(Y));
1753: PetscFunctionReturn(PETSC_SUCCESS);
1754: }
1756: /*@
1757: VecCreateLocalVector - Creates a vector object suitable for use with `VecGetLocalVector()` and friends. You must call `VecDestroy()` when the
1758: vector is no longer needed.
1760: Not Collective.
1762: Input Parameter:
1763: . v - The vector for which the local vector is desired.
1765: Output Parameter:
1766: . w - Upon exit this contains the local vector.
1768: Level: beginner
1770: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1771: @*/
1772: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1773: {
1774: PetscMPIInt size;
1776: PetscFunctionBegin;
1778: PetscAssertPointer(w, 2);
1779: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1780: if (size == 1) PetscCall(VecDuplicate(v, w));
1781: else if (v->ops->createlocalvector) PetscUseTypeMethod(v, createlocalvector, w);
1782: else {
1783: VecType type;
1784: PetscInt n;
1786: PetscCall(VecCreate(PETSC_COMM_SELF, w));
1787: PetscCall(VecGetLocalSize(v, &n));
1788: PetscCall(VecSetSizes(*w, n, n));
1789: PetscCall(VecGetBlockSize(v, &n));
1790: PetscCall(VecSetBlockSize(*w, n));
1791: PetscCall(VecGetType(v, &type));
1792: PetscCall(VecSetType(*w, type));
1793: }
1794: PetscFunctionReturn(PETSC_SUCCESS);
1795: }
1797: /*@
1798: VecGetLocalVectorRead - Maps the local portion of a vector into a
1799: vector.
1801: Not Collective.
1803: Input Parameter:
1804: . v - The vector for which the local vector is desired.
1806: Output Parameter:
1807: . w - Upon exit this contains the local vector.
1809: Level: beginner
1811: Notes:
1812: You must call `VecRestoreLocalVectorRead()` when the local
1813: vector is no longer needed.
1815: This function is similar to `VecGetArrayRead()` which maps the local
1816: portion into a raw pointer. `VecGetLocalVectorRead()` is usually
1817: almost as efficient as `VecGetArrayRead()` but in certain circumstances
1818: `VecGetLocalVectorRead()` can be much more efficient than
1819: `VecGetArrayRead()`. This is because the construction of a contiguous
1820: array representing the vector data required by `VecGetArrayRead()` can
1821: be an expensive operation for certain vector types. For example, for
1822: GPU vectors `VecGetArrayRead()` requires that the data between device
1823: and host is synchronized.
1825: Unlike `VecGetLocalVector()`, this routine is not collective and
1826: preserves cached information.
1828: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1829: @*/
1830: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1831: {
1832: PetscFunctionBegin;
1835: VecCheckSameLocalSize(v, 1, w, 2);
1836: if (v->ops->getlocalvectorread) {
1837: PetscUseTypeMethod(v, getlocalvectorread, w);
1838: } else {
1839: PetscScalar *a;
1841: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1842: PetscCall(VecPlaceArray(w, a));
1843: }
1844: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1845: PetscCall(VecLockReadPush(v));
1846: PetscCall(VecLockReadPush(w));
1847: PetscFunctionReturn(PETSC_SUCCESS);
1848: }
1850: /*@
1851: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1852: previously mapped into a vector using `VecGetLocalVectorRead()`.
1854: Not Collective.
1856: Input Parameters:
1857: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
1858: - w - The vector into which the local portion of `v` was mapped.
1860: Level: beginner
1862: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1863: @*/
1864: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1865: {
1866: PetscFunctionBegin;
1869: if (v->ops->restorelocalvectorread) {
1870: PetscUseTypeMethod(v, restorelocalvectorread, w);
1871: } else {
1872: const PetscScalar *a;
1874: PetscCall(VecGetArrayRead(w, &a));
1875: PetscCall(VecRestoreArrayRead(v, &a));
1876: PetscCall(VecResetArray(w));
1877: }
1878: PetscCall(VecLockReadPop(v));
1879: PetscCall(VecLockReadPop(w));
1880: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1881: PetscFunctionReturn(PETSC_SUCCESS);
1882: }
1884: /*@
1885: VecGetLocalVector - Maps the local portion of a vector into a
1886: vector.
1888: Collective
1890: Input Parameter:
1891: . v - The vector for which the local vector is desired.
1893: Output Parameter:
1894: . w - Upon exit this contains the local vector.
1896: Level: beginner
1898: Notes:
1899: You must call `VecRestoreLocalVector()` when the local
1900: vector is no longer needed.
1902: This function is similar to `VecGetArray()` which maps the local
1903: portion into a raw pointer. `VecGetLocalVector()` is usually about as
1904: efficient as `VecGetArray()` but in certain circumstances
1905: `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1906: This is because the construction of a contiguous array representing
1907: the vector data required by `VecGetArray()` can be an expensive
1908: operation for certain vector types. For example, for GPU vectors
1909: `VecGetArray()` requires that the data between device and host is
1910: synchronized.
1912: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1913: @*/
1914: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1915: {
1916: PetscFunctionBegin;
1919: VecCheckSameLocalSize(v, 1, w, 2);
1920: if (v->ops->getlocalvector) {
1921: PetscUseTypeMethod(v, getlocalvector, w);
1922: } else {
1923: PetscScalar *a;
1925: PetscCall(VecGetArray(v, &a));
1926: PetscCall(VecPlaceArray(w, a));
1927: }
1928: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1929: PetscFunctionReturn(PETSC_SUCCESS);
1930: }
1932: /*@
1933: VecRestoreLocalVector - Unmaps the local portion of a vector
1934: previously mapped into a vector using `VecGetLocalVector()`.
1936: Logically Collective.
1938: Input Parameters:
1939: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
1940: - w - The vector into which the local portion of `v` was mapped.
1942: Level: beginner
1944: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1945: @*/
1946: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1947: {
1948: PetscFunctionBegin;
1951: if (v->ops->restorelocalvector) {
1952: PetscUseTypeMethod(v, restorelocalvector, w);
1953: } else {
1954: PetscScalar *a;
1955: PetscCall(VecGetArray(w, &a));
1956: PetscCall(VecRestoreArray(v, &a));
1957: PetscCall(VecResetArray(w));
1958: }
1959: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1960: PetscCall(PetscObjectStateIncrease((PetscObject)v));
1961: PetscFunctionReturn(PETSC_SUCCESS);
1962: }
1964: /*@C
1965: VecGetArray - Returns a pointer to a contiguous array that contains this
1966: MPI processes's portion of the vector data
1968: Logically Collective
1970: Input Parameter:
1971: . x - the vector
1973: Output Parameter:
1974: . a - location to put pointer to the array
1976: Level: beginner
1978: Notes:
1979: For the standard PETSc vectors, `VecGetArray()` returns a pointer to the local data array and
1980: does not use any copies. If the underlying vector data is not stored in a contiguous array
1981: this routine will copy the data to a contiguous array and return a pointer to that. You MUST
1982: call `VecRestoreArray()` when you no longer need access to the array.
1984: Fortran Notes:
1985: `VecGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayF90()`
1987: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
1988: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
1989: @*/
1990: PetscErrorCode VecGetArray(Vec x, PetscScalar **a)
1991: {
1992: PetscFunctionBegin;
1994: PetscCall(VecSetErrorIfLocked(x, 1));
1995: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
1996: PetscUseTypeMethod(x, getarray, a);
1997: } else if (x->petscnative) { /* VECSTANDARD */
1998: *a = *((PetscScalar **)x->data);
1999: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
2000: PetscFunctionReturn(PETSC_SUCCESS);
2001: }
2003: /*@C
2004: VecRestoreArray - Restores a vector after `VecGetArray()` has been called and the array is no longer needed
2006: Logically Collective
2008: Input Parameters:
2009: + x - the vector
2010: - a - location of pointer to array obtained from `VecGetArray()`
2012: Level: beginner
2014: Fortran Notes:
2015: `VecRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayF90()`
2017: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2018: `VecGetArrayPair()`, `VecRestoreArrayPair()`
2019: @*/
2020: PetscErrorCode VecRestoreArray(Vec x, PetscScalar **a)
2021: {
2022: PetscFunctionBegin;
2024: if (a) PetscAssertPointer(a, 2);
2025: if (x->ops->restorearray) {
2026: PetscUseTypeMethod(x, restorearray, a);
2027: } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
2028: if (a) *a = NULL;
2029: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2030: PetscFunctionReturn(PETSC_SUCCESS);
2031: }
2032: /*@C
2033: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
2035: Not Collective
2037: Input Parameter:
2038: . x - the vector
2040: Output Parameter:
2041: . a - the array
2043: Level: beginner
2045: Notes:
2046: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2048: Unlike `VecGetArray()`, preserves cached information like vector norms.
2050: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
2051: implementations may require a copy, but such implementations should cache the contiguous representation so that
2052: only one copy is performed when this routine is called multiple times in sequence.
2054: Fortran Notes:
2055: `VecGetArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayReadF90()`
2057: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2058: @*/
2059: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar **a)
2060: {
2061: PetscFunctionBegin;
2063: PetscAssertPointer(a, 2);
2064: if (x->ops->getarrayread) {
2065: PetscUseTypeMethod(x, getarrayread, a);
2066: } else if (x->ops->getarray) {
2067: PetscObjectState state;
2069: /* VECNEST, VECCUDA, VECKOKKOS etc */
2070: // x->ops->getarray may bump the object state, but since we know this is a read-only get
2071: // we can just undo that
2072: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2073: PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
2074: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2075: } else if (x->petscnative) {
2076: /* VECSTANDARD */
2077: *a = *((PetscScalar **)x->data);
2078: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
2079: PetscFunctionReturn(PETSC_SUCCESS);
2080: }
2082: /*@C
2083: VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`
2085: Not Collective
2087: Input Parameters:
2088: + x - the vector
2089: - a - the array
2091: Level: beginner
2093: Fortran Notes:
2094: `VecRestoreArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayReadF90()`
2096: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2097: @*/
2098: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar **a)
2099: {
2100: PetscFunctionBegin;
2102: if (a) PetscAssertPointer(a, 2);
2103: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
2104: /* nothing */
2105: } else if (x->ops->restorearrayread) { /* VECNEST */
2106: PetscUseTypeMethod(x, restorearrayread, a);
2107: } else { /* No one? */
2108: PetscObjectState state;
2110: // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2111: // we can just undo that
2112: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2113: PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2114: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2115: }
2116: if (a) *a = NULL;
2117: PetscFunctionReturn(PETSC_SUCCESS);
2118: }
2120: /*@C
2121: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
2122: MPI processes's portion of the vector data.
2124: Logically Collective
2126: Input Parameter:
2127: . x - the vector
2129: Output Parameter:
2130: . a - location to put pointer to the array
2132: Level: intermediate
2134: Note:
2135: The values in this array are NOT valid, the caller of this routine is responsible for putting
2136: values into the array; any values it does not set will be invalid.
2138: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2140: For vectors associated with GPUs, the host and device vectors are not synchronized before
2141: giving access. If you need correct values in the array use `VecGetArray()`
2143: Fortran Notes:
2144: `VecGetArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayWriteF90()`
2146: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteF90()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
2147: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
2148: @*/
2149: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar **a)
2150: {
2151: PetscFunctionBegin;
2153: PetscAssertPointer(a, 2);
2154: PetscCall(VecSetErrorIfLocked(x, 1));
2155: if (x->ops->getarraywrite) {
2156: PetscUseTypeMethod(x, getarraywrite, a);
2157: } else {
2158: PetscCall(VecGetArray(x, a));
2159: }
2160: PetscFunctionReturn(PETSC_SUCCESS);
2161: }
2163: /*@C
2164: VecRestoreArrayWrite - Restores a vector after `VecGetArrayWrite()` has been called.
2166: Logically Collective
2168: Input Parameters:
2169: + x - the vector
2170: - a - location of pointer to array obtained from `VecGetArray()`
2172: Level: beginner
2174: Fortran Notes:
2175: `VecRestoreArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayWriteF90()`
2177: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteF90()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2178: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2179: @*/
2180: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar **a)
2181: {
2182: PetscFunctionBegin;
2184: if (a) PetscAssertPointer(a, 2);
2185: if (x->ops->restorearraywrite) {
2186: PetscUseTypeMethod(x, restorearraywrite, a);
2187: } else if (x->ops->restorearray) {
2188: PetscUseTypeMethod(x, restorearray, a);
2189: }
2190: if (a) *a = NULL;
2191: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2192: PetscFunctionReturn(PETSC_SUCCESS);
2193: }
2195: /*@C
2196: VecGetArrays - Returns a pointer to the arrays in a set of vectors
2197: that were created by a call to `VecDuplicateVecs()`.
2199: Logically Collective; No Fortran Support
2201: Input Parameters:
2202: + x - the vectors
2203: - n - the number of vectors
2205: Output Parameter:
2206: . a - location to put pointer to the array
2208: Level: intermediate
2210: Note:
2211: You MUST call `VecRestoreArrays()` when you no longer need access to the arrays.
2213: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2214: @*/
2215: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2216: {
2217: PetscInt i;
2218: PetscScalar **q;
2220: PetscFunctionBegin;
2221: PetscAssertPointer(x, 1);
2223: PetscAssertPointer(a, 3);
2224: PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2225: PetscCall(PetscMalloc1(n, &q));
2226: for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2227: *a = q;
2228: PetscFunctionReturn(PETSC_SUCCESS);
2229: }
2231: /*@C
2232: VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2233: has been called.
2235: Logically Collective; No Fortran Support
2237: Input Parameters:
2238: + x - the vector
2239: . n - the number of vectors
2240: - a - location of pointer to arrays obtained from `VecGetArrays()`
2242: Notes:
2243: For regular PETSc vectors this routine does not involve any copies. For
2244: any special vectors that do not store local vector data in a contiguous
2245: array, this routine will copy the data back into the underlying
2246: vector data structure from the arrays obtained with `VecGetArrays()`.
2248: Level: intermediate
2250: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2251: @*/
2252: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2253: {
2254: PetscInt i;
2255: PetscScalar **q = *a;
2257: PetscFunctionBegin;
2258: PetscAssertPointer(x, 1);
2260: PetscAssertPointer(a, 3);
2262: for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2263: PetscCall(PetscFree(q));
2264: PetscFunctionReturn(PETSC_SUCCESS);
2265: }
2267: /*@C
2268: VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g.,
2269: `VECCUDA`), the returned pointer will be a device pointer to the device memory that contains
2270: this MPI processes's portion of the vector data.
2272: Logically Collective; No Fortran Support
2274: Input Parameter:
2275: . x - the vector
2277: Output Parameters:
2278: + a - location to put pointer to the array
2279: - mtype - memory type of the array
2281: Level: beginner
2283: Note:
2284: Device data is guaranteed to have the latest value. Otherwise, when this is a host vector
2285: (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host
2286: pointer.
2288: For `VECKOKKOS`, if Kokkos is configured without device (e.g., use serial or openmp), per
2289: this function, the vector works like `VECSEQ`/`VECMPI`; otherwise, it works like `VECCUDA` or
2290: `VECHIP` etc.
2292: Use `VecRestoreArrayAndMemType()` when the array access is no longer needed.
2294: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`,
2295: `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2296: @*/
2297: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2298: {
2299: PetscFunctionBegin;
2302: PetscAssertPointer(a, 2);
2303: if (mtype) PetscAssertPointer(mtype, 3);
2304: PetscCall(VecSetErrorIfLocked(x, 1));
2305: if (x->ops->getarrayandmemtype) {
2306: /* VECCUDA, VECKOKKOS etc */
2307: PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2308: } else {
2309: /* VECSTANDARD, VECNEST, VECVIENNACL */
2310: PetscCall(VecGetArray(x, a));
2311: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2312: }
2313: PetscFunctionReturn(PETSC_SUCCESS);
2314: }
2316: /*@C
2317: VecRestoreArrayAndMemType - Restores a vector after `VecGetArrayAndMemType()` has been called.
2319: Logically Collective; No Fortran Support
2321: Input Parameters:
2322: + x - the vector
2323: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`
2325: Level: beginner
2327: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`,
2328: `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2329: @*/
2330: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar **a)
2331: {
2332: PetscFunctionBegin;
2335: if (a) PetscAssertPointer(a, 2);
2336: if (x->ops->restorearrayandmemtype) {
2337: /* VECCUDA, VECKOKKOS etc */
2338: PetscUseTypeMethod(x, restorearrayandmemtype, a);
2339: } else {
2340: /* VECNEST, VECVIENNACL */
2341: PetscCall(VecRestoreArray(x, a));
2342: } /* VECSTANDARD does nothing */
2343: if (a) *a = NULL;
2344: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2345: PetscFunctionReturn(PETSC_SUCCESS);
2346: }
2348: /*@C
2349: VecGetArrayReadAndMemType - Like `VecGetArrayRead()`, but if the input vector is a device vector, it will return a read-only device pointer.
2350: The returned pointer is guaranteed to point to up-to-date data. For host vectors, it functions as `VecGetArrayRead()`.
2352: Not Collective; No Fortran Support
2354: Input Parameter:
2355: . x - the vector
2357: Output Parameters:
2358: + a - the array
2359: - mtype - memory type of the array
2361: Level: beginner
2363: Notes:
2364: The array must be returned using a matching call to `VecRestoreArrayReadAndMemType()`.
2366: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2367: @*/
2368: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar **a, PetscMemType *mtype)
2369: {
2370: PetscFunctionBegin;
2373: PetscAssertPointer(a, 2);
2374: if (mtype) PetscAssertPointer(mtype, 3);
2375: if (x->ops->getarrayreadandmemtype) {
2376: /* VECCUDA/VECHIP though they are also petscnative */
2377: PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2378: } else if (x->ops->getarrayandmemtype) {
2379: /* VECKOKKOS */
2380: PetscObjectState state;
2382: // see VecGetArrayRead() for why
2383: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2384: PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2385: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2386: } else {
2387: PetscCall(VecGetArrayRead(x, a));
2388: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2389: }
2390: PetscFunctionReturn(PETSC_SUCCESS);
2391: }
2393: /*@C
2394: VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`
2396: Not Collective; No Fortran Support
2398: Input Parameters:
2399: + x - the vector
2400: - a - the array
2402: Level: beginner
2404: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2405: @*/
2406: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar **a)
2407: {
2408: PetscFunctionBegin;
2411: if (a) PetscAssertPointer(a, 2);
2412: if (x->ops->restorearrayreadandmemtype) {
2413: /* VECCUDA/VECHIP */
2414: PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2415: } else if (!x->petscnative) {
2416: /* VECNEST */
2417: PetscCall(VecRestoreArrayRead(x, a));
2418: }
2419: if (a) *a = NULL;
2420: PetscFunctionReturn(PETSC_SUCCESS);
2421: }
2423: /*@C
2424: VecGetArrayWriteAndMemType - Like `VecGetArrayWrite()`, but if this is a device vector it will always return
2425: a device pointer to the device memory that contains this processor's portion of the vector data.
2427: Logically Collective; No Fortran Support
2429: Input Parameter:
2430: . x - the vector
2432: Output Parameters:
2433: + a - the array
2434: - mtype - memory type of the array
2436: Level: beginner
2438: Note:
2439: The array must be returned using a matching call to `VecRestoreArrayWriteAndMemType()`, where it will label the device memory as most recent.
2441: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2442: @*/
2443: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2444: {
2445: PetscFunctionBegin;
2448: PetscCall(VecSetErrorIfLocked(x, 1));
2449: PetscAssertPointer(a, 2);
2450: if (mtype) PetscAssertPointer(mtype, 3);
2451: if (x->ops->getarraywriteandmemtype) {
2452: /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2453: PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2454: } else if (x->ops->getarrayandmemtype) {
2455: PetscCall(VecGetArrayAndMemType(x, a, mtype));
2456: } else {
2457: /* VECNEST, VECVIENNACL */
2458: PetscCall(VecGetArrayWrite(x, a));
2459: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2460: }
2461: PetscFunctionReturn(PETSC_SUCCESS);
2462: }
2464: /*@C
2465: VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`
2467: Logically Collective; No Fortran Support
2469: Input Parameters:
2470: + x - the vector
2471: - a - the array
2473: Level: beginner
2475: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2476: @*/
2477: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar **a)
2478: {
2479: PetscFunctionBegin;
2482: PetscCall(VecSetErrorIfLocked(x, 1));
2483: if (a) PetscAssertPointer(a, 2);
2484: if (x->ops->restorearraywriteandmemtype) {
2485: /* VECCUDA/VECHIP */
2486: PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2487: PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2488: } else if (x->ops->restorearrayandmemtype) {
2489: PetscCall(VecRestoreArrayAndMemType(x, a));
2490: } else {
2491: PetscCall(VecRestoreArray(x, a));
2492: }
2493: if (a) *a = NULL;
2494: PetscFunctionReturn(PETSC_SUCCESS);
2495: }
2497: /*@
2498: VecPlaceArray - Allows one to replace the array in a vector with an
2499: array provided by the user. This is useful to avoid copying an array
2500: into a vector.
2502: Logically Collective; No Fortran Support
2504: Input Parameters:
2505: + vec - the vector
2506: - array - the array
2508: Level: developer
2510: Notes:
2511: Use `VecReplaceArray()` instead to permanently replace the array
2513: You can return to the original array with a call to `VecResetArray()`. `vec` does not take
2514: ownership of `array` in any way.
2516: The user must free `array` themselves but be careful not to
2517: do so before the vector has either been destroyed, had its original array restored with
2518: `VecResetArray()` or permanently replaced with `VecReplaceArray()`.
2520: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2521: @*/
2522: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2523: {
2524: PetscFunctionBegin;
2527: if (array) PetscAssertPointer(array, 2);
2528: PetscUseTypeMethod(vec, placearray, array);
2529: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2530: PetscFunctionReturn(PETSC_SUCCESS);
2531: }
2533: /*@C
2534: VecReplaceArray - Allows one to replace the array in a vector with an
2535: array provided by the user. This is useful to avoid copying an array
2536: into a vector.
2538: Logically Collective; No Fortran Support
2540: Input Parameters:
2541: + vec - the vector
2542: - array - the array
2544: Level: developer
2546: Notes:
2547: This permanently replaces the array and frees the memory associated
2548: with the old array. Use `VecPlaceArray()` to temporarily replace the array.
2550: The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
2551: freed by the user. It will be freed when the vector is destroyed.
2553: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2554: @*/
2555: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2556: {
2557: PetscFunctionBegin;
2560: PetscUseTypeMethod(vec, replacearray, array);
2561: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2562: PetscFunctionReturn(PETSC_SUCCESS);
2563: }
2565: /*MC
2566: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
2567: and makes them accessible via a Fortran pointer.
2569: Synopsis:
2570: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
2572: Collective
2574: Input Parameters:
2575: + x - a vector to mimic
2576: - n - the number of vectors to obtain
2578: Output Parameters:
2579: + y - Fortran pointer to the array of vectors
2580: - ierr - error code
2582: Example of Usage:
2583: .vb
2584: #include <petsc/finclude/petscvec.h>
2585: use petscvec
2587: Vec x
2588: Vec, pointer :: y(:)
2589: ....
2590: call VecDuplicateVecsF90(x,2,y,ierr)
2591: call VecSet(y(2),alpha,ierr)
2592: call VecSet(y(2),alpha,ierr)
2593: ....
2594: call VecDestroyVecsF90(2,y,ierr)
2595: .ve
2597: Level: beginner
2599: Note:
2600: Use `VecDestroyVecsF90()` to free the space.
2602: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecsF90()`, `VecDuplicateVecs()`
2603: M*/
2605: /*MC
2606: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2607: `VecGetArrayF90()`.
2609: Synopsis:
2610: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2612: Logically Collective
2614: Input Parameters:
2615: + x - vector
2616: - xx_v - the Fortran pointer to the array
2618: Output Parameter:
2619: . ierr - error code
2621: Example of Usage:
2622: .vb
2623: #include <petsc/finclude/petscvec.h>
2624: use petscvec
2626: PetscScalar, pointer :: xx_v(:)
2627: ....
2628: call VecGetArrayF90(x,xx_v,ierr)
2629: xx_v(3) = a
2630: call VecRestoreArrayF90(x,xx_v,ierr)
2631: .ve
2633: Level: beginner
2635: .seealso: [](ch_vectors), `Vec`, `VecGetArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrayReadF90()`
2636: M*/
2638: /*MC
2639: VecDestroyVecsF90 - Frees a block of vectors obtained with `VecDuplicateVecsF90()`.
2641: Synopsis:
2642: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
2644: Collective
2646: Input Parameters:
2647: + n - the number of vectors previously obtained
2648: - x - pointer to array of vector pointers
2650: Output Parameter:
2651: . ierr - error code
2653: Level: beginner
2655: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecs()`, `VecDuplicateVecsF90()`
2656: M*/
2658: /*MC
2659: VecGetArrayF90 - Accesses a vector array from Fortran. For default PETSc
2660: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2661: this routine is implementation dependent. You MUST call `VecRestoreArrayF90()`
2662: when you no longer need access to the array.
2664: Synopsis:
2665: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2667: Logically Collective
2669: Input Parameter:
2670: . x - vector
2672: Output Parameters:
2673: + xx_v - the Fortran pointer to the array
2674: - ierr - error code
2676: Example of Usage:
2677: .vb
2678: #include <petsc/finclude/petscvec.h>
2679: use petscvec
2681: PetscScalar, pointer :: xx_v(:)
2682: ....
2683: call VecGetArrayF90(x,xx_v,ierr)
2684: xx_v(3) = a
2685: call VecRestoreArrayF90(x,xx_v,ierr)
2686: .ve
2688: Level: beginner
2690: Note:
2691: If you ONLY intend to read entries from the array and not change any entries you should use `VecGetArrayReadF90()`.
2693: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayReadF90()`
2694: M*/
2696: /*MC
2697: VecGetArrayReadF90 - Accesses a read only array from Fortran. For default PETSc
2698: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2699: this routine is implementation dependent. You MUST call `VecRestoreArrayReadF90()`
2700: when you no longer need access to the array.
2702: Synopsis:
2703: VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2705: Logically Collective
2707: Input Parameter:
2708: . x - vector
2710: Output Parameters:
2711: + xx_v - the Fortran pointer to the array
2712: - ierr - error code
2714: Example of Usage:
2715: .vb
2716: #include <petsc/finclude/petscvec.h>
2717: use petscvec
2719: PetscScalar, pointer :: xx_v(:)
2720: ....
2721: call VecGetArrayReadF90(x,xx_v,ierr)
2722: a = xx_v(3)
2723: call VecRestoreArrayReadF90(x,xx_v,ierr)
2724: .ve
2726: Level: beginner
2728: Note:
2729: If you intend to write entries into the array you must use `VecGetArrayF90()`.
2731: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecGetArrayF90()`
2732: M*/
2734: /*MC
2735: VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2736: `VecGetArrayReadF90()`.
2738: Synopsis:
2739: VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2741: Logically Collective
2743: Input Parameters:
2744: + x - vector
2745: - xx_v - the Fortran pointer to the array
2747: Output Parameter:
2748: . ierr - error code
2750: Example of Usage:
2751: .vb
2752: #include <petsc/finclude/petscvec.h>
2753: use petscvec
2755: PetscScalar, pointer :: xx_v(:)
2756: ....
2757: call VecGetArrayReadF90(x,xx_v,ierr)
2758: a = xx_v(3)
2759: call VecRestoreArrayReadF90(x,xx_v,ierr)
2760: .ve
2762: Level: beginner
2764: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecRestoreArrayF90()`
2765: M*/
2767: /*@C
2768: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2769: processor's portion of the vector data. You MUST call `VecRestoreArray2d()`
2770: when you no longer need access to the array.
2772: Logically Collective
2774: Input Parameters:
2775: + x - the vector
2776: . m - first dimension of two dimensional array
2777: . n - second dimension of two dimensional array
2778: . mstart - first index you will use in first coordinate direction (often 0)
2779: - nstart - first index in the second coordinate direction (often 0)
2781: Output Parameter:
2782: . a - location to put pointer to the array
2784: Level: developer
2786: Notes:
2787: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2788: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2789: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2790: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2792: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2794: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2795: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2796: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2797: @*/
2798: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2799: {
2800: PetscInt i, N;
2801: PetscScalar *aa;
2803: PetscFunctionBegin;
2805: PetscAssertPointer(a, 6);
2807: PetscCall(VecGetLocalSize(x, &N));
2808: 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);
2809: PetscCall(VecGetArray(x, &aa));
2811: PetscCall(PetscMalloc1(m, a));
2812: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2813: *a -= mstart;
2814: PetscFunctionReturn(PETSC_SUCCESS);
2815: }
2817: /*@C
2818: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2819: processor's portion of the vector data. You MUST call `VecRestoreArray2dWrite()`
2820: when you no longer need access to the array.
2822: Logically Collective
2824: Input Parameters:
2825: + x - the vector
2826: . m - first dimension of two dimensional array
2827: . n - second dimension of two dimensional array
2828: . mstart - first index you will use in first coordinate direction (often 0)
2829: - nstart - first index in the second coordinate direction (often 0)
2831: Output Parameter:
2832: . a - location to put pointer to the array
2834: Level: developer
2836: Notes:
2837: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2838: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2839: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2840: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2842: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2844: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2845: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2846: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2847: @*/
2848: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2849: {
2850: PetscInt i, N;
2851: PetscScalar *aa;
2853: PetscFunctionBegin;
2855: PetscAssertPointer(a, 6);
2857: PetscCall(VecGetLocalSize(x, &N));
2858: 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);
2859: PetscCall(VecGetArrayWrite(x, &aa));
2861: PetscCall(PetscMalloc1(m, a));
2862: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2863: *a -= mstart;
2864: PetscFunctionReturn(PETSC_SUCCESS);
2865: }
2867: /*@C
2868: VecRestoreArray2d - Restores a vector after `VecGetArray2d()` has been called.
2870: Logically Collective
2872: Input Parameters:
2873: + x - the vector
2874: . m - first dimension of two dimensional array
2875: . n - second dimension of the two dimensional array
2876: . mstart - first index you will use in first coordinate direction (often 0)
2877: . nstart - first index in the second coordinate direction (often 0)
2878: - a - location of pointer to array obtained from `VecGetArray2d()`
2880: Level: developer
2882: Notes:
2883: For regular PETSc vectors this routine does not involve any copies. For
2884: any special vectors that do not store local vector data in a contiguous
2885: array, this routine will copy the data back into the underlying
2886: vector data structure from the array obtained with `VecGetArray()`.
2888: This routine actually zeros out the `a` pointer.
2890: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2891: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2892: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2893: @*/
2894: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2895: {
2896: void *dummy;
2898: PetscFunctionBegin;
2900: PetscAssertPointer(a, 6);
2902: dummy = (void *)(*a + mstart);
2903: PetscCall(PetscFree(dummy));
2904: PetscCall(VecRestoreArray(x, NULL));
2905: *a = NULL;
2906: PetscFunctionReturn(PETSC_SUCCESS);
2907: }
2909: /*@C
2910: VecRestoreArray2dWrite - Restores a vector after `VecGetArray2dWrite()` has been called.
2912: Logically Collective
2914: Input Parameters:
2915: + x - the vector
2916: . m - first dimension of two dimensional array
2917: . n - second dimension of the two dimensional array
2918: . mstart - first index you will use in first coordinate direction (often 0)
2919: . nstart - first index in the second coordinate direction (often 0)
2920: - a - location of pointer to array obtained from `VecGetArray2d()`
2922: Level: developer
2924: Notes:
2925: For regular PETSc vectors this routine does not involve any copies. For
2926: any special vectors that do not store local vector data in a contiguous
2927: array, this routine will copy the data back into the underlying
2928: vector data structure from the array obtained with `VecGetArray()`.
2930: This routine actually zeros out the `a` pointer.
2932: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2933: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2934: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2935: @*/
2936: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2937: {
2938: void *dummy;
2940: PetscFunctionBegin;
2942: PetscAssertPointer(a, 6);
2944: dummy = (void *)(*a + mstart);
2945: PetscCall(PetscFree(dummy));
2946: PetscCall(VecRestoreArrayWrite(x, NULL));
2947: PetscFunctionReturn(PETSC_SUCCESS);
2948: }
2950: /*@C
2951: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2952: processor's portion of the vector data. You MUST call `VecRestoreArray1d()`
2953: when you no longer need access to the array.
2955: Logically Collective
2957: Input Parameters:
2958: + x - the vector
2959: . m - first dimension of two dimensional array
2960: - mstart - first index you will use in first coordinate direction (often 0)
2962: Output Parameter:
2963: . a - location to put pointer to the array
2965: Level: developer
2967: Notes:
2968: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2969: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2970: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2972: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2974: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2975: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2976: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2977: @*/
2978: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2979: {
2980: PetscInt N;
2982: PetscFunctionBegin;
2984: PetscAssertPointer(a, 4);
2986: PetscCall(VecGetLocalSize(x, &N));
2987: 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);
2988: PetscCall(VecGetArray(x, a));
2989: *a -= mstart;
2990: PetscFunctionReturn(PETSC_SUCCESS);
2991: }
2993: /*@C
2994: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
2995: processor's portion of the vector data. You MUST call `VecRestoreArray1dWrite()`
2996: when you no longer need access to the array.
2998: Logically Collective
3000: Input Parameters:
3001: + x - the vector
3002: . m - first dimension of two dimensional array
3003: - mstart - first index you will use in first coordinate direction (often 0)
3005: Output Parameter:
3006: . a - location to put pointer to the array
3008: Level: developer
3010: Notes:
3011: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3012: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3013: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3015: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3017: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3018: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3019: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3020: @*/
3021: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3022: {
3023: PetscInt N;
3025: PetscFunctionBegin;
3027: PetscAssertPointer(a, 4);
3029: PetscCall(VecGetLocalSize(x, &N));
3030: 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);
3031: PetscCall(VecGetArrayWrite(x, a));
3032: *a -= mstart;
3033: PetscFunctionReturn(PETSC_SUCCESS);
3034: }
3036: /*@C
3037: VecRestoreArray1d - Restores a vector after `VecGetArray1d()` has been called.
3039: Logically Collective
3041: Input Parameters:
3042: + x - the vector
3043: . m - first dimension of two dimensional array
3044: . mstart - first index you will use in first coordinate direction (often 0)
3045: - a - location of pointer to array obtained from `VecGetArray1d()`
3047: Level: developer
3049: Notes:
3050: For regular PETSc vectors this routine does not involve any copies. For
3051: any special vectors that do not store local vector data in a contiguous
3052: array, this routine will copy the data back into the underlying
3053: vector data structure from the array obtained with `VecGetArray1d()`.
3055: This routine actually zeros out the `a` pointer.
3057: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3058: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3059: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3060: @*/
3061: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3062: {
3063: PetscFunctionBegin;
3066: PetscCall(VecRestoreArray(x, NULL));
3067: *a = NULL;
3068: PetscFunctionReturn(PETSC_SUCCESS);
3069: }
3071: /*@C
3072: VecRestoreArray1dWrite - Restores a vector after `VecGetArray1dWrite()` has been called.
3074: Logically Collective
3076: Input Parameters:
3077: + x - the vector
3078: . m - first dimension of two dimensional array
3079: . mstart - first index you will use in first coordinate direction (often 0)
3080: - a - location of pointer to array obtained from `VecGetArray1d()`
3082: Level: developer
3084: Notes:
3085: For regular PETSc vectors this routine does not involve any copies. For
3086: any special vectors that do not store local vector data in a contiguous
3087: array, this routine will copy the data back into the underlying
3088: vector data structure from the array obtained with `VecGetArray1d()`.
3090: This routine actually zeros out the `a` pointer.
3092: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3093: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3094: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3095: @*/
3096: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3097: {
3098: PetscFunctionBegin;
3101: PetscCall(VecRestoreArrayWrite(x, NULL));
3102: *a = NULL;
3103: PetscFunctionReturn(PETSC_SUCCESS);
3104: }
3106: /*@C
3107: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
3108: processor's portion of the vector data. You MUST call `VecRestoreArray3d()`
3109: when you no longer need access to the array.
3111: Logically Collective
3113: Input Parameters:
3114: + x - the vector
3115: . m - first dimension of three dimensional array
3116: . n - second dimension of three dimensional array
3117: . p - third dimension of three dimensional array
3118: . mstart - first index you will use in first coordinate direction (often 0)
3119: . nstart - first index in the second coordinate direction (often 0)
3120: - pstart - first index in the third coordinate direction (often 0)
3122: Output Parameter:
3123: . a - location to put pointer to the array
3125: Level: developer
3127: Notes:
3128: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3129: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3130: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3131: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3133: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3135: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3136: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
3137: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3138: @*/
3139: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3140: {
3141: PetscInt i, N, j;
3142: PetscScalar *aa, **b;
3144: PetscFunctionBegin;
3146: PetscAssertPointer(a, 8);
3148: PetscCall(VecGetLocalSize(x, &N));
3149: 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);
3150: PetscCall(VecGetArray(x, &aa));
3152: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3153: b = (PetscScalar **)((*a) + m);
3154: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3155: for (i = 0; i < m; i++)
3156: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset(aa, i * n * p + j * p - pstart);
3157: *a -= mstart;
3158: PetscFunctionReturn(PETSC_SUCCESS);
3159: }
3161: /*@C
3162: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
3163: processor's portion of the vector data. You MUST call `VecRestoreArray3dWrite()`
3164: when you no longer need access to the array.
3166: Logically Collective
3168: Input Parameters:
3169: + x - the vector
3170: . m - first dimension of three dimensional array
3171: . n - second dimension of three dimensional array
3172: . p - third dimension of three dimensional array
3173: . mstart - first index you will use in first coordinate direction (often 0)
3174: . nstart - first index in the second coordinate direction (often 0)
3175: - pstart - first index in the third coordinate direction (often 0)
3177: Output Parameter:
3178: . a - location to put pointer to the array
3180: Level: developer
3182: Notes:
3183: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3184: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3185: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3186: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3188: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3190: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3191: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3192: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3193: @*/
3194: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3195: {
3196: PetscInt i, N, j;
3197: PetscScalar *aa, **b;
3199: PetscFunctionBegin;
3201: PetscAssertPointer(a, 8);
3203: PetscCall(VecGetLocalSize(x, &N));
3204: 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);
3205: PetscCall(VecGetArrayWrite(x, &aa));
3207: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3208: b = (PetscScalar **)((*a) + m);
3209: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3210: for (i = 0; i < m; i++)
3211: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3213: *a -= mstart;
3214: PetscFunctionReturn(PETSC_SUCCESS);
3215: }
3217: /*@C
3218: VecRestoreArray3d - Restores a vector after `VecGetArray3d()` has been called.
3220: Logically Collective
3222: Input Parameters:
3223: + x - the vector
3224: . m - first dimension of three dimensional array
3225: . n - second dimension of the three dimensional array
3226: . p - third dimension of the three dimensional array
3227: . mstart - first index you will use in first coordinate direction (often 0)
3228: . nstart - first index in the second coordinate direction (often 0)
3229: . pstart - first index in the third coordinate direction (often 0)
3230: - a - location of pointer to array obtained from VecGetArray3d()
3232: Level: developer
3234: Notes:
3235: For regular PETSc vectors this routine does not involve any copies. For
3236: any special vectors that do not store local vector data in a contiguous
3237: array, this routine will copy the data back into the underlying
3238: vector data structure from the array obtained with `VecGetArray()`.
3240: This routine actually zeros out the `a` pointer.
3242: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3243: `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3244: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3245: @*/
3246: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3247: {
3248: void *dummy;
3250: PetscFunctionBegin;
3252: PetscAssertPointer(a, 8);
3254: dummy = (void *)(*a + mstart);
3255: PetscCall(PetscFree(dummy));
3256: PetscCall(VecRestoreArray(x, NULL));
3257: *a = NULL;
3258: PetscFunctionReturn(PETSC_SUCCESS);
3259: }
3261: /*@C
3262: VecRestoreArray3dWrite - Restores a vector after `VecGetArray3dWrite()` has been called.
3264: Logically Collective
3266: Input Parameters:
3267: + x - the vector
3268: . m - first dimension of three dimensional array
3269: . n - second dimension of the three dimensional array
3270: . p - third dimension of the three dimensional array
3271: . mstart - first index you will use in first coordinate direction (often 0)
3272: . nstart - first index in the second coordinate direction (often 0)
3273: . pstart - first index in the third coordinate direction (often 0)
3274: - a - location of pointer to array obtained from VecGetArray3d()
3276: Level: developer
3278: Notes:
3279: For regular PETSc vectors this routine does not involve any copies. For
3280: any special vectors that do not store local vector data in a contiguous
3281: array, this routine will copy the data back into the underlying
3282: vector data structure from the array obtained with `VecGetArray()`.
3284: This routine actually zeros out the `a` pointer.
3286: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3287: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3288: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3289: @*/
3290: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3291: {
3292: void *dummy;
3294: PetscFunctionBegin;
3296: PetscAssertPointer(a, 8);
3298: dummy = (void *)(*a + mstart);
3299: PetscCall(PetscFree(dummy));
3300: PetscCall(VecRestoreArrayWrite(x, NULL));
3301: *a = NULL;
3302: PetscFunctionReturn(PETSC_SUCCESS);
3303: }
3305: /*@C
3306: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3307: processor's portion of the vector data. You MUST call `VecRestoreArray4d()`
3308: when you no longer need access to the array.
3310: Logically Collective
3312: Input Parameters:
3313: + x - the vector
3314: . m - first dimension of four dimensional array
3315: . n - second dimension of four dimensional array
3316: . p - third dimension of four dimensional array
3317: . q - fourth dimension of four dimensional array
3318: . mstart - first index you will use in first coordinate direction (often 0)
3319: . nstart - first index in the second coordinate direction (often 0)
3320: . pstart - first index in the third coordinate direction (often 0)
3321: - qstart - first index in the fourth coordinate direction (often 0)
3323: Output Parameter:
3324: . a - location to put pointer to the array
3326: Level: developer
3328: Notes:
3329: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3330: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3331: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3332: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3334: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3336: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3337: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3338: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3339: @*/
3340: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3341: {
3342: PetscInt i, N, j, k;
3343: PetscScalar *aa, ***b, **c;
3345: PetscFunctionBegin;
3347: PetscAssertPointer(a, 10);
3349: PetscCall(VecGetLocalSize(x, &N));
3350: 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);
3351: PetscCall(VecGetArray(x, &aa));
3353: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3354: b = (PetscScalar ***)((*a) + m);
3355: c = (PetscScalar **)(b + m * n);
3356: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3357: for (i = 0; i < m; i++)
3358: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3359: for (i = 0; i < m; i++)
3360: for (j = 0; j < n; j++)
3361: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3362: *a -= mstart;
3363: PetscFunctionReturn(PETSC_SUCCESS);
3364: }
3366: /*@C
3367: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3368: processor's portion of the vector data. You MUST call `VecRestoreArray4dWrite()`
3369: when you no longer need access to the array.
3371: Logically Collective
3373: Input Parameters:
3374: + x - the vector
3375: . m - first dimension of four dimensional array
3376: . n - second dimension of four dimensional array
3377: . p - third dimension of four dimensional array
3378: . q - fourth dimension of four dimensional array
3379: . mstart - first index you will use in first coordinate direction (often 0)
3380: . nstart - first index in the second coordinate direction (often 0)
3381: . pstart - first index in the third coordinate direction (often 0)
3382: - qstart - first index in the fourth coordinate direction (often 0)
3384: Output Parameter:
3385: . a - location to put pointer to the array
3387: Level: developer
3389: Notes:
3390: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3391: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3392: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3393: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3395: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3397: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3398: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3399: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3400: @*/
3401: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3402: {
3403: PetscInt i, N, j, k;
3404: PetscScalar *aa, ***b, **c;
3406: PetscFunctionBegin;
3408: PetscAssertPointer(a, 10);
3410: PetscCall(VecGetLocalSize(x, &N));
3411: 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);
3412: PetscCall(VecGetArrayWrite(x, &aa));
3414: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3415: b = (PetscScalar ***)((*a) + m);
3416: c = (PetscScalar **)(b + m * n);
3417: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3418: for (i = 0; i < m; i++)
3419: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3420: for (i = 0; i < m; i++)
3421: for (j = 0; j < n; j++)
3422: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3423: *a -= mstart;
3424: PetscFunctionReturn(PETSC_SUCCESS);
3425: }
3427: /*@C
3428: VecRestoreArray4d - Restores a vector after `VecGetArray4d()` has been called.
3430: Logically Collective
3432: Input Parameters:
3433: + x - the vector
3434: . m - first dimension of four dimensional array
3435: . n - second dimension of the four dimensional array
3436: . p - third dimension of the four dimensional array
3437: . q - fourth dimension of the four dimensional array
3438: . mstart - first index you will use in first coordinate direction (often 0)
3439: . nstart - first index in the second coordinate direction (often 0)
3440: . pstart - first index in the third coordinate direction (often 0)
3441: . qstart - first index in the fourth coordinate direction (often 0)
3442: - a - location of pointer to array obtained from VecGetArray4d()
3444: Level: developer
3446: Notes:
3447: For regular PETSc vectors this routine does not involve any copies. For
3448: any special vectors that do not store local vector data in a contiguous
3449: array, this routine will copy the data back into the underlying
3450: vector data structure from the array obtained with `VecGetArray()`.
3452: This routine actually zeros out the `a` pointer.
3454: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3455: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3456: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecGet`
3457: @*/
3458: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3459: {
3460: void *dummy;
3462: PetscFunctionBegin;
3464: PetscAssertPointer(a, 10);
3466: dummy = (void *)(*a + mstart);
3467: PetscCall(PetscFree(dummy));
3468: PetscCall(VecRestoreArray(x, NULL));
3469: *a = NULL;
3470: PetscFunctionReturn(PETSC_SUCCESS);
3471: }
3473: /*@C
3474: VecRestoreArray4dWrite - Restores a vector after `VecGetArray4dWrite()` has been called.
3476: Logically Collective
3478: Input Parameters:
3479: + x - the vector
3480: . m - first dimension of four dimensional array
3481: . n - second dimension of the four dimensional array
3482: . p - third dimension of the four dimensional array
3483: . q - fourth dimension of the four dimensional array
3484: . mstart - first index you will use in first coordinate direction (often 0)
3485: . nstart - first index in the second coordinate direction (often 0)
3486: . pstart - first index in the third coordinate direction (often 0)
3487: . qstart - first index in the fourth coordinate direction (often 0)
3488: - a - location of pointer to array obtained from `VecGetArray4d()`
3490: Level: developer
3492: Notes:
3493: For regular PETSc vectors this routine does not involve any copies. For
3494: any special vectors that do not store local vector data in a contiguous
3495: array, this routine will copy the data back into the underlying
3496: vector data structure from the array obtained with `VecGetArray()`.
3498: This routine actually zeros out the `a` pointer.
3500: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3501: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3502: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3503: @*/
3504: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3505: {
3506: void *dummy;
3508: PetscFunctionBegin;
3510: PetscAssertPointer(a, 10);
3512: dummy = (void *)(*a + mstart);
3513: PetscCall(PetscFree(dummy));
3514: PetscCall(VecRestoreArrayWrite(x, NULL));
3515: *a = NULL;
3516: PetscFunctionReturn(PETSC_SUCCESS);
3517: }
3519: /*@C
3520: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3521: processor's portion of the vector data. You MUST call `VecRestoreArray2dRead()`
3522: when you no longer need access to the array.
3524: Logically Collective
3526: Input Parameters:
3527: + x - the vector
3528: . m - first dimension of two dimensional array
3529: . n - second dimension of two dimensional array
3530: . mstart - first index you will use in first coordinate direction (often 0)
3531: - nstart - first index in the second coordinate direction (often 0)
3533: Output Parameter:
3534: . a - location to put pointer to the array
3536: Level: developer
3538: Notes:
3539: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
3540: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3541: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3542: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
3544: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3546: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3547: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3548: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3549: @*/
3550: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3551: {
3552: PetscInt i, N;
3553: const PetscScalar *aa;
3555: PetscFunctionBegin;
3557: PetscAssertPointer(a, 6);
3559: PetscCall(VecGetLocalSize(x, &N));
3560: 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);
3561: PetscCall(VecGetArrayRead(x, &aa));
3563: PetscCall(PetscMalloc1(m, a));
3564: for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3565: *a -= mstart;
3566: PetscFunctionReturn(PETSC_SUCCESS);
3567: }
3569: /*@C
3570: VecRestoreArray2dRead - Restores a vector after `VecGetArray2dRead()` has been called.
3572: Logically Collective
3574: Input Parameters:
3575: + x - the vector
3576: . m - first dimension of two dimensional array
3577: . n - second dimension of the two dimensional array
3578: . mstart - first index you will use in first coordinate direction (often 0)
3579: . nstart - first index in the second coordinate direction (often 0)
3580: - a - location of pointer to array obtained from VecGetArray2d()
3582: Level: developer
3584: Notes:
3585: For regular PETSc vectors this routine does not involve any copies. For
3586: any special vectors that do not store local vector data in a contiguous
3587: array, this routine will copy the data back into the underlying
3588: vector data structure from the array obtained with `VecGetArray()`.
3590: This routine actually zeros out the `a` pointer.
3592: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3593: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3594: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3595: @*/
3596: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3597: {
3598: void *dummy;
3600: PetscFunctionBegin;
3602: PetscAssertPointer(a, 6);
3604: dummy = (void *)(*a + mstart);
3605: PetscCall(PetscFree(dummy));
3606: PetscCall(VecRestoreArrayRead(x, NULL));
3607: *a = NULL;
3608: PetscFunctionReturn(PETSC_SUCCESS);
3609: }
3611: /*@C
3612: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3613: processor's portion of the vector data. You MUST call `VecRestoreArray1dRead()`
3614: when you no longer need access to the array.
3616: Logically Collective
3618: Input Parameters:
3619: + x - the vector
3620: . m - first dimension of two dimensional array
3621: - mstart - first index you will use in first coordinate direction (often 0)
3623: Output Parameter:
3624: . a - location to put pointer to the array
3626: Level: developer
3628: Notes:
3629: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3630: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3631: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3633: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3635: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3636: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3637: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3638: @*/
3639: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3640: {
3641: PetscInt N;
3643: PetscFunctionBegin;
3645: PetscAssertPointer(a, 4);
3647: PetscCall(VecGetLocalSize(x, &N));
3648: 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);
3649: PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3650: *a -= mstart;
3651: PetscFunctionReturn(PETSC_SUCCESS);
3652: }
3654: /*@C
3655: VecRestoreArray1dRead - Restores a vector after `VecGetArray1dRead()` has been called.
3657: Logically Collective
3659: Input Parameters:
3660: + x - the vector
3661: . m - first dimension of two dimensional array
3662: . mstart - first index you will use in first coordinate direction (often 0)
3663: - a - location of pointer to array obtained from `VecGetArray1dRead()`
3665: Level: developer
3667: Notes:
3668: For regular PETSc vectors this routine does not involve any copies. For
3669: any special vectors that do not store local vector data in a contiguous
3670: array, this routine will copy the data back into the underlying
3671: vector data structure from the array obtained with `VecGetArray1dRead()`.
3673: This routine actually zeros out the `a` pointer.
3675: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3676: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3677: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3678: @*/
3679: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3680: {
3681: PetscFunctionBegin;
3684: PetscCall(VecRestoreArrayRead(x, NULL));
3685: *a = NULL;
3686: PetscFunctionReturn(PETSC_SUCCESS);
3687: }
3689: /*@C
3690: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3691: processor's portion of the vector data. You MUST call `VecRestoreArray3dRead()`
3692: when you no longer need access to the array.
3694: Logically Collective
3696: Input Parameters:
3697: + x - the vector
3698: . m - first dimension of three dimensional array
3699: . n - second dimension of three dimensional array
3700: . p - third dimension of three dimensional array
3701: . mstart - first index you will use in first coordinate direction (often 0)
3702: . nstart - first index in the second coordinate direction (often 0)
3703: - pstart - first index in the third coordinate direction (often 0)
3705: Output Parameter:
3706: . a - location to put pointer to the array
3708: Level: developer
3710: Notes:
3711: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3712: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3713: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3714: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3dRead()`.
3716: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3718: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3719: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3720: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3721: @*/
3722: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3723: {
3724: PetscInt i, N, j;
3725: const PetscScalar *aa;
3726: PetscScalar **b;
3728: PetscFunctionBegin;
3730: PetscAssertPointer(a, 8);
3732: PetscCall(VecGetLocalSize(x, &N));
3733: 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);
3734: PetscCall(VecGetArrayRead(x, &aa));
3736: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3737: b = (PetscScalar **)((*a) + m);
3738: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3739: for (i = 0; i < m; i++)
3740: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset((PetscScalar *)aa, i * n * p + j * p - pstart);
3741: *a -= mstart;
3742: PetscFunctionReturn(PETSC_SUCCESS);
3743: }
3745: /*@C
3746: VecRestoreArray3dRead - Restores a vector after `VecGetArray3dRead()` has been called.
3748: Logically Collective
3750: Input Parameters:
3751: + x - the vector
3752: . m - first dimension of three dimensional array
3753: . n - second dimension of the three dimensional array
3754: . p - third dimension of the three dimensional array
3755: . mstart - first index you will use in first coordinate direction (often 0)
3756: . nstart - first index in the second coordinate direction (often 0)
3757: . pstart - first index in the third coordinate direction (often 0)
3758: - a - location of pointer to array obtained from `VecGetArray3dRead()`
3760: Level: developer
3762: Notes:
3763: For regular PETSc vectors this routine does not involve any copies. For
3764: any special vectors that do not store local vector data in a contiguous
3765: array, this routine will copy the data back into the underlying
3766: vector data structure from the array obtained with `VecGetArray()`.
3768: This routine actually zeros out the `a` pointer.
3770: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3771: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3772: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3773: @*/
3774: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3775: {
3776: void *dummy;
3778: PetscFunctionBegin;
3780: PetscAssertPointer(a, 8);
3782: dummy = (void *)(*a + mstart);
3783: PetscCall(PetscFree(dummy));
3784: PetscCall(VecRestoreArrayRead(x, NULL));
3785: *a = NULL;
3786: PetscFunctionReturn(PETSC_SUCCESS);
3787: }
3789: /*@C
3790: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3791: processor's portion of the vector data. You MUST call `VecRestoreArray4dRead()`
3792: when you no longer need access to the array.
3794: Logically Collective
3796: Input Parameters:
3797: + x - the vector
3798: . m - first dimension of four dimensional array
3799: . n - second dimension of four dimensional array
3800: . p - third dimension of four dimensional array
3801: . q - fourth dimension of four dimensional array
3802: . mstart - first index you will use in first coordinate direction (often 0)
3803: . nstart - first index in the second coordinate direction (often 0)
3804: . pstart - first index in the third coordinate direction (often 0)
3805: - qstart - first index in the fourth coordinate direction (often 0)
3807: Output Parameter:
3808: . a - location to put pointer to the array
3810: Level: beginner
3812: Notes:
3813: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3814: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3815: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3816: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3818: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3820: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3821: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3822: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3823: @*/
3824: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3825: {
3826: PetscInt i, N, j, k;
3827: const PetscScalar *aa;
3828: PetscScalar ***b, **c;
3830: PetscFunctionBegin;
3832: PetscAssertPointer(a, 10);
3834: PetscCall(VecGetLocalSize(x, &N));
3835: 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);
3836: PetscCall(VecGetArrayRead(x, &aa));
3838: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3839: b = (PetscScalar ***)((*a) + m);
3840: c = (PetscScalar **)(b + m * n);
3841: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3842: for (i = 0; i < m; i++)
3843: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3844: for (i = 0; i < m; i++)
3845: for (j = 0; j < n; j++)
3846: 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;
3847: *a -= mstart;
3848: PetscFunctionReturn(PETSC_SUCCESS);
3849: }
3851: /*@C
3852: VecRestoreArray4dRead - Restores a vector after `VecGetArray4d()` has been called.
3854: Logically Collective
3856: Input Parameters:
3857: + x - the vector
3858: . m - first dimension of four dimensional array
3859: . n - second dimension of the four dimensional array
3860: . p - third dimension of the four dimensional array
3861: . q - fourth dimension of the four dimensional array
3862: . mstart - first index you will use in first coordinate direction (often 0)
3863: . nstart - first index in the second coordinate direction (often 0)
3864: . pstart - first index in the third coordinate direction (often 0)
3865: . qstart - first index in the fourth coordinate direction (often 0)
3866: - a - location of pointer to array obtained from `VecGetArray4dRead()`
3868: Level: beginner
3870: Notes:
3871: For regular PETSc vectors this routine does not involve any copies. For
3872: any special vectors that do not store local vector data in a contiguous
3873: array, this routine will copy the data back into the underlying
3874: vector data structure from the array obtained with `VecGetArray()`.
3876: This routine actually zeros out the `a` pointer.
3878: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3879: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3880: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3881: @*/
3882: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3883: {
3884: void *dummy;
3886: PetscFunctionBegin;
3888: PetscAssertPointer(a, 10);
3890: dummy = (void *)(*a + mstart);
3891: PetscCall(PetscFree(dummy));
3892: PetscCall(VecRestoreArrayRead(x, NULL));
3893: *a = NULL;
3894: PetscFunctionReturn(PETSC_SUCCESS);
3895: }
3897: #if defined(PETSC_USE_DEBUG)
3899: /*@
3900: VecLockGet - Gets the current lock status of a vector
3902: Logically Collective
3904: Input Parameter:
3905: . x - the vector
3907: Output Parameter:
3908: . state - greater than zero indicates the vector is locked for read; less than zero indicates the vector is
3909: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3911: Level: advanced
3913: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3914: @*/
3915: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3916: {
3917: PetscFunctionBegin;
3919: *state = x->lock;
3920: PetscFunctionReturn(PETSC_SUCCESS);
3921: }
3923: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3924: {
3925: PetscFunctionBegin;
3927: PetscAssertPointer(file, 2);
3928: PetscAssertPointer(func, 3);
3929: PetscAssertPointer(line, 4);
3930: #if !PetscDefined(HAVE_THREADSAFETY)
3931: {
3932: const int index = x->lockstack.currentsize - 1;
3934: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupted vec lock stack, have negative index %d", index);
3935: *file = x->lockstack.file[index];
3936: *func = x->lockstack.function[index];
3937: *line = x->lockstack.line[index];
3938: }
3939: #else
3940: *file = NULL;
3941: *func = NULL;
3942: *line = 0;
3943: #endif
3944: PetscFunctionReturn(PETSC_SUCCESS);
3945: }
3947: /*@
3948: VecLockReadPush - Pushes a read-only lock on a vector to prevent it from being written to
3950: Logically Collective
3952: Input Parameter:
3953: . x - the vector
3955: Level: intermediate
3957: Notes:
3958: If this is set then calls to `VecGetArray()` or `VecSetValues()` or any other routines that change the vectors values will generate an error.
3960: The call can be nested, i.e., called multiple times on the same vector, but each `VecLockReadPush()` has to have one matching
3961: `VecLockReadPop()`, which removes the latest read-only lock.
3963: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3964: @*/
3965: PetscErrorCode VecLockReadPush(Vec x)
3966: {
3967: PetscFunctionBegin;
3969: 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");
3970: #if !PetscDefined(HAVE_THREADSAFETY)
3971: {
3972: const char *file, *func;
3973: int index, line;
3975: if ((index = petscstack.currentsize - 2) == -1) {
3976: // vec was locked "outside" of petsc, either in user-land or main. the error message will
3977: // now show this function as the culprit, but it will include the stacktrace
3978: file = "unknown user-file";
3979: func = "unknown_user_function";
3980: line = 0;
3981: } else {
3982: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected petscstack, have negative index %d", index);
3983: file = petscstack.file[index];
3984: func = petscstack.function[index];
3985: line = petscstack.line[index];
3986: }
3987: PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
3988: }
3989: #endif
3990: PetscFunctionReturn(PETSC_SUCCESS);
3991: }
3993: /*@
3994: VecLockReadPop - Pops a read-only lock from a vector
3996: Logically Collective
3998: Input Parameter:
3999: . x - the vector
4001: Level: intermediate
4003: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
4004: @*/
4005: PetscErrorCode VecLockReadPop(Vec x)
4006: {
4007: PetscFunctionBegin;
4009: PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
4010: #if !PetscDefined(HAVE_THREADSAFETY)
4011: {
4012: const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];
4014: PetscStackPop_Private(x->lockstack, previous);
4015: }
4016: #endif
4017: PetscFunctionReturn(PETSC_SUCCESS);
4018: }
4020: /*@
4021: VecLockWriteSet - Lock or unlock a vector for exclusive read/write access
4023: Logically Collective
4025: Input Parameters:
4026: + x - the vector
4027: - flg - `PETSC_TRUE` to lock the vector for exclusive read/write access; `PETSC_FALSE` to unlock it.
4029: Level: intermediate
4031: Notes:
4032: The function is useful in split-phase computations, which usually have a begin phase and an end phase.
4033: One can call `VecLockWriteSet`(x,`PETSC_TRUE`) in the begin phase to lock a vector for exclusive
4034: access, and call `VecLockWriteSet`(x,`PETSC_FALSE`) in the end phase to unlock the vector from exclusive
4035: access. In this way, one is ensured no other operations can access the vector in between. The code may like
4037: .vb
4038: VecGetArray(x,&xdata); // begin phase
4039: VecLockWriteSet(v,PETSC_TRUE);
4041: Other operations, which can not access x anymore (they can access xdata, of course)
4043: VecRestoreArray(x,&vdata); // end phase
4044: VecLockWriteSet(v,PETSC_FALSE);
4045: .ve
4047: The call can not be nested on the same vector, in other words, one can not call `VecLockWriteSet`(x,`PETSC_TRUE`)
4048: again before calling `VecLockWriteSet`(v,`PETSC_FALSE`).
4050: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
4051: @*/
4052: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
4053: {
4054: PetscFunctionBegin;
4056: if (flg) {
4057: 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");
4058: 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");
4059: x->lock = -1;
4060: } else {
4061: 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");
4062: x->lock = 0;
4063: }
4064: PetscFunctionReturn(PETSC_SUCCESS);
4065: }
4067: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
4068: /*@
4069: VecLockPush - Pushes a read-only lock on a vector to prevent it from being written to
4071: Level: deprecated
4073: .seealso: [](ch_vectors), `Vec`, `VecLockReadPush()`
4074: @*/
4075: PetscErrorCode VecLockPush(Vec x)
4076: {
4077: PetscFunctionBegin;
4078: PetscCall(VecLockReadPush(x));
4079: PetscFunctionReturn(PETSC_SUCCESS);
4080: }
4082: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
4083: /*@
4084: VecLockPop - Pops a read-only lock from a vector
4086: Level: deprecated
4088: .seealso: [](ch_vectors), `Vec`, `VecLockReadPop()`
4089: @*/
4090: PetscErrorCode VecLockPop(Vec x)
4091: {
4092: PetscFunctionBegin;
4093: PetscCall(VecLockReadPop(x));
4094: PetscFunctionReturn(PETSC_SUCCESS);
4095: }
4097: #endif