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: PetscCallMPI(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: PetscCallMPI(MPIU_Allreduce(b1, b2, 2, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)x)));
219: PetscCheck((PetscIsNanReal(b2[0]) && PetscIsNanReal(b2[1])) || (-b2[0] == b2[1]), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Difference in cached %s norms: local %g", NormTypes[type], (double)*val);
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: /*@
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: Fortran Note:
911: If any of `ix` and `y` are scalars pass them using, for example,
912: .vb
913: VecSetValues(mat, one, [ix], [y], INSERT_VALUES)
914: .ve
916: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
917: `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
918: @*/
919: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
920: {
921: PetscFunctionBeginHot;
923: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
924: PetscAssertPointer(ix, 3);
925: PetscAssertPointer(y, 4);
928: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
929: PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
930: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
931: PetscCall(PetscObjectStateIncrease((PetscObject)x));
932: PetscFunctionReturn(PETSC_SUCCESS);
933: }
935: /*@
936: VecGetValues - Gets values from certain locations of a vector. Currently
937: can only get values on the same processor on which they are owned
939: Not Collective
941: Input Parameters:
942: + x - vector to get values from
943: . ni - number of elements to get
944: - ix - indices where to get them from (in global 1d numbering)
946: Output Parameter:
947: . y - array of values, must be passed in with a length of `ni`
949: Level: beginner
951: Notes:
952: The user provides the allocated array y; it is NOT allocated in this routine
954: `VecGetValues()` gets y[i] = x[ix[i]], for i=0,...,ni-1.
956: `VecAssemblyBegin()` and `VecAssemblyEnd()` MUST be called before calling this if `VecSetValues()` or related routine has been called
958: VecGetValues() uses 0-based indices in Fortran as well as in C.
960: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
961: negative indices may be passed in ix. These rows are
962: simply ignored.
964: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
965: @*/
966: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
967: {
968: PetscFunctionBegin;
970: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
971: PetscAssertPointer(ix, 3);
972: PetscAssertPointer(y, 4);
974: VecCheckAssembled(x);
975: PetscUseTypeMethod(x, getvalues, ni, ix, y);
976: PetscFunctionReturn(PETSC_SUCCESS);
977: }
979: /*@
980: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
982: Not Collective
984: Input Parameters:
985: + x - vector to insert in
986: . ni - number of blocks to add
987: . ix - indices where to add in block count, rather than element count
988: . y - array of values
989: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries
991: Level: intermediate
993: Notes:
994: `VecSetValuesBlocked()` sets x[bs*ix[i]+j] = y[bs*i+j],
995: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
997: Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
998: options cannot be mixed without intervening calls to the assembly
999: routines.
1001: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1002: MUST be called after all calls to `VecSetValuesBlocked()` have been completed.
1004: `VecSetValuesBlocked()` uses 0-based indices in Fortran as well as in C.
1006: Negative indices may be passed in ix, these rows are
1007: simply ignored. This allows easily inserting element load matrices
1008: with homogeneous Dirichlet boundary conditions that you don't want represented
1009: in the vector.
1011: Fortran Note:
1012: If any of `ix` and `y` are scalars pass them using, for example,
1013: .vb
1014: VecSetValuesBlocked(mat, one, [ix], [y], INSERT_VALUES)
1015: .ve
1017: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
1018: `VecSetValues()`
1019: @*/
1020: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1021: {
1022: PetscFunctionBeginHot;
1024: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1025: PetscAssertPointer(ix, 3);
1026: PetscAssertPointer(y, 4);
1029: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1030: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1031: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1032: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1033: PetscFunctionReturn(PETSC_SUCCESS);
1034: }
1036: /*@
1037: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1038: using a local ordering of the nodes.
1040: Not Collective
1042: Input Parameters:
1043: + x - vector to insert in
1044: . ni - number of elements to add
1045: . ix - indices where to add
1046: . y - array of values
1047: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1049: Level: intermediate
1051: Notes:
1052: `VecSetValuesLocal()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
1054: Calls to `VecSetValuesLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1055: options cannot be mixed without intervening calls to the assembly
1056: routines.
1058: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1059: MUST be called after all calls to `VecSetValuesLocal()` have been completed.
1061: `VecSetValuesLocal()` uses 0-based indices in Fortran as well as in C.
1063: Fortran Note:
1064: If any of `ix` and `y` are scalars pass them using, for example,
1065: .vb
1066: VecSetValuesLocal(mat, one, [ix], [y], INSERT_VALUES)
1067: .ve
1069: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1070: `VecSetValuesBlockedLocal()`
1071: @*/
1072: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1073: {
1074: PetscInt lixp[128], *lix = lixp;
1076: PetscFunctionBeginHot;
1078: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1079: PetscAssertPointer(ix, 3);
1080: PetscAssertPointer(y, 4);
1083: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1084: if (!x->ops->setvalueslocal) {
1085: if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
1086: if (x->map->mapping) {
1087: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1088: PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1089: PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1090: if (ni > 128) PetscCall(PetscFree(lix));
1091: } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1092: } else PetscUseTypeMethod(x, setvalueslocal, ni, ix, y, iora);
1093: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1094: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1095: PetscFunctionReturn(PETSC_SUCCESS);
1096: }
1098: /*@
1099: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1100: using a local ordering of the nodes.
1102: Not Collective
1104: Input Parameters:
1105: + x - vector to insert in
1106: . ni - number of blocks to add
1107: . ix - indices where to add in block count, not element count
1108: . y - array of values
1109: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1111: Level: intermediate
1113: Notes:
1114: `VecSetValuesBlockedLocal()` sets x[bs*ix[i]+j] = y[bs*i+j],
1115: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with `VecSetBlockSize()`.
1117: Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1118: options cannot be mixed without intervening calls to the assembly
1119: routines.
1121: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1122: MUST be called after all calls to `VecSetValuesBlockedLocal()` have been completed.
1124: `VecSetValuesBlockedLocal()` uses 0-based indices in Fortran as well as in C.
1126: Fortran Note:
1127: If any of `ix` and `y` are scalars pass them using, for example,
1128: .vb
1129: VecSetValuesBlockedLocal(mat, one, [ix], [y], INSERT_VALUES)
1130: .ve
1132: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1133: `VecSetLocalToGlobalMapping()`
1134: @*/
1135: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1136: {
1137: PetscInt lixp[128], *lix = lixp;
1139: PetscFunctionBeginHot;
1141: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1142: PetscAssertPointer(ix, 3);
1143: PetscAssertPointer(y, 4);
1145: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1146: if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
1147: if (x->map->mapping) {
1148: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1149: PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1150: PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1151: if (ni > 128) PetscCall(PetscFree(lix));
1152: } else {
1153: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1154: }
1155: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1156: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1157: PetscFunctionReturn(PETSC_SUCCESS);
1158: }
1160: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1161: {
1162: PetscFunctionBegin;
1165: VecCheckAssembled(x);
1167: if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1168: PetscAssertPointer(y, 3);
1169: for (PetscInt i = 0; i < nv; ++i) {
1172: PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1173: VecCheckSameSize(x, 1, y[i], 3);
1174: VecCheckAssembled(y[i]);
1175: PetscCall(VecLockReadPush(y[i]));
1176: }
1177: PetscAssertPointer(result, 4);
1180: PetscCall(VecLockReadPush(x));
1181: PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1182: PetscCall((*mxdot)(x, nv, y, result));
1183: PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1184: PetscCall(VecLockReadPop(x));
1185: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1186: PetscFunctionReturn(PETSC_SUCCESS);
1187: }
1189: /*@
1190: VecMTDot - Computes indefinite vector multiple dot products.
1191: That is, it does NOT use the complex conjugate.
1193: Collective
1195: Input Parameters:
1196: + x - one vector
1197: . nv - number of vectors
1198: - y - array of vectors. Note that vectors are pointers
1200: Output Parameter:
1201: . val - array of the dot products
1203: Level: intermediate
1205: Notes for Users of Complex Numbers:
1206: For complex vectors, `VecMTDot()` computes the indefinite form
1207: $ val = (x,y) = y^T x,
1208: where y^T denotes the transpose of y.
1210: Use `VecMDot()` for the inner product
1211: $ val = (x,y) = y^H x,
1212: where y^H denotes the conjugate transpose of y.
1214: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1215: @*/
1216: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1217: {
1218: PetscFunctionBegin;
1220: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1221: PetscFunctionReturn(PETSC_SUCCESS);
1222: }
1224: /*@
1225: VecMDot - Computes multiple vector dot products.
1227: Collective
1229: Input Parameters:
1230: + x - one vector
1231: . nv - number of vectors
1232: - y - array of vectors.
1234: Output Parameter:
1235: . val - array of the dot products (does not allocate the array)
1237: Level: intermediate
1239: Notes for Users of Complex Numbers:
1240: For complex vectors, `VecMDot()` computes
1241: $ val = (x,y) = y^H x,
1242: where y^H denotes the conjugate transpose of y.
1244: Use `VecMTDot()` for the indefinite form
1245: $ val = (x,y) = y^T x,
1246: where y^T denotes the transpose of y.
1248: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`
1249: @*/
1250: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1251: {
1252: PetscFunctionBegin;
1254: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1255: PetscFunctionReturn(PETSC_SUCCESS);
1256: }
1258: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
1259: {
1260: PetscFunctionBegin;
1262: VecCheckAssembled(y);
1264: PetscCall(VecSetErrorIfLocked(y, 1));
1265: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1266: if (nv) {
1267: PetscInt zeros = 0;
1269: PetscAssertPointer(alpha, 3);
1270: PetscAssertPointer(x, 4);
1271: for (PetscInt i = 0; i < nv; ++i) {
1275: PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1276: VecCheckSameSize(y, 1, x[i], 4);
1277: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1278: VecCheckAssembled(x[i]);
1279: PetscCall(VecLockReadPush(x[i]));
1280: zeros += alpha[i] == (PetscScalar)0.0;
1281: }
1283: if (zeros < nv) {
1284: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1285: VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1286: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1287: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1288: }
1290: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1291: }
1292: PetscFunctionReturn(PETSC_SUCCESS);
1293: }
1295: /*@
1296: VecMAXPY - Computes `y = y + sum alpha[i] x[i]`
1298: Logically Collective
1300: Input Parameters:
1301: + nv - number of scalars and x-vectors
1302: . alpha - array of scalars
1303: . y - one vector
1304: - x - array of vectors
1306: Level: intermediate
1308: Note:
1309: `y` cannot be any of the `x` vectors
1311: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`,`VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1312: @*/
1313: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1314: {
1315: PetscFunctionBegin;
1316: PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1317: PetscFunctionReturn(PETSC_SUCCESS);
1318: }
1320: /*@
1321: VecMAXPBY - Computes `y = beta y + sum alpha[i] x[i]`
1323: Logically Collective
1325: Input Parameters:
1326: + nv - number of scalars and x-vectors
1327: . alpha - array of scalars
1328: . beta - scalar
1329: . y - one vector
1330: - x - array of vectors
1332: Level: intermediate
1334: Note:
1335: `y` cannot be any of the `x` vectors.
1337: Developer Notes:
1338: This is a convenience routine, but implementations might be able to optimize it, for example, when `beta` is zero.
1340: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1341: @*/
1342: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
1343: {
1344: PetscFunctionBegin;
1346: VecCheckAssembled(y);
1348: PetscCall(VecSetErrorIfLocked(y, 1));
1349: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1352: if (y->ops->maxpby) {
1353: PetscInt zeros = 0;
1355: if (nv) {
1356: PetscAssertPointer(alpha, 3);
1357: PetscAssertPointer(x, 5);
1358: }
1360: for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
1364: PetscCheckSameTypeAndComm(y, 1, x[i], 5);
1365: VecCheckSameSize(y, 1, x[i], 5);
1366: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1367: VecCheckAssembled(x[i]);
1368: PetscCall(VecLockReadPush(x[i]));
1369: zeros += alpha[i] == (PetscScalar)0.0;
1370: }
1372: if (zeros < nv) { // has nonzero alpha
1373: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1374: PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1375: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1376: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1377: } else {
1378: PetscCall(VecScale(y, beta));
1379: }
1381: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1382: } else { // no maxpby
1383: if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1384: else PetscCall(VecScale(y, beta));
1385: PetscCall(VecMAXPY(y, nv, alpha, x));
1386: }
1387: PetscFunctionReturn(PETSC_SUCCESS);
1388: }
1390: /*@
1391: VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
1392: in the order they appear in the array. The concatenated vector resides on the same
1393: communicator and is the same type as the source vectors.
1395: Collective
1397: Input Parameters:
1398: + nx - number of vectors to be concatenated
1399: - X - array containing the vectors to be concatenated in the order of concatenation
1401: Output Parameters:
1402: + Y - concatenated vector
1403: - x_is - array of index sets corresponding to the concatenated components of `Y` (pass `NULL` if not needed)
1405: Level: advanced
1407: Notes:
1408: Concatenation is similar to the functionality of a `VECNEST` object; they both represent combination of
1409: different vector spaces. However, concatenated vectors do not store any information about their
1410: sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1411: manipulation of data in the concatenated vector that corresponds to the original components at creation.
1413: This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1414: has to operate on combined vector spaces and cannot utilize `VECNEST` objects due to incompatibility with
1415: bound projections.
1417: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1418: @*/
1419: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1420: {
1421: MPI_Comm comm;
1422: VecType vec_type;
1423: Vec Ytmp, Xtmp;
1424: IS *is_tmp;
1425: PetscInt i, shift = 0, Xnl, Xng, Xbegin;
1427: PetscFunctionBegin;
1431: PetscAssertPointer(Y, 3);
1433: if ((*X)->ops->concatenate) {
1434: /* use the dedicated concatenation function if available */
1435: PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1436: } else {
1437: /* loop over vectors and start creating IS */
1438: comm = PetscObjectComm((PetscObject)*X);
1439: PetscCall(VecGetType(*X, &vec_type));
1440: PetscCall(PetscMalloc1(nx, &is_tmp));
1441: for (i = 0; i < nx; i++) {
1442: PetscCall(VecGetSize(X[i], &Xng));
1443: PetscCall(VecGetLocalSize(X[i], &Xnl));
1444: PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1445: PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1446: shift += Xng;
1447: }
1448: /* create the concatenated vector */
1449: PetscCall(VecCreate(comm, &Ytmp));
1450: PetscCall(VecSetType(Ytmp, vec_type));
1451: PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1452: PetscCall(VecSetUp(Ytmp));
1453: /* copy data from X array to Y and return */
1454: for (i = 0; i < nx; i++) {
1455: PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1456: PetscCall(VecCopy(X[i], Xtmp));
1457: PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1458: }
1459: *Y = Ytmp;
1460: if (x_is) {
1461: *x_is = is_tmp;
1462: } else {
1463: for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1464: PetscCall(PetscFree(is_tmp));
1465: }
1466: }
1467: PetscFunctionReturn(PETSC_SUCCESS);
1468: }
1470: /* A helper function for VecGetSubVector to check if we can implement it with no-copy (i.e. the subvector shares
1471: memory with the original vector), and the block size of the subvector.
1473: Input Parameters:
1474: + X - the original vector
1475: - is - the index set of the subvector
1477: Output Parameters:
1478: + contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
1479: . start - start of contiguous block, as an offset from the start of the ownership range of the original vector
1480: - blocksize - the block size of the subvector
1482: */
1483: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1484: {
1485: PetscInt gstart, gend, lstart;
1486: PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1487: PetscInt n, N, ibs, vbs, bs = -1;
1489: PetscFunctionBegin;
1490: PetscCall(ISGetLocalSize(is, &n));
1491: PetscCall(ISGetSize(is, &N));
1492: PetscCall(ISGetBlockSize(is, &ibs));
1493: PetscCall(VecGetBlockSize(X, &vbs));
1494: PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1495: PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1496: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1497: if (ibs > 1) {
1498: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1499: bs = ibs;
1500: } else {
1501: if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1502: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1503: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1504: }
1506: *contig = red[0];
1507: *start = lstart;
1508: *blocksize = bs;
1509: PetscFunctionReturn(PETSC_SUCCESS);
1510: }
1512: /* A helper function for VecGetSubVector, to be used when we have to build a standalone subvector through VecScatter
1514: Input Parameters:
1515: + X - the original vector
1516: . is - the index set of the subvector
1517: - bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()
1519: Output Parameter:
1520: . Z - the subvector, which will compose the VecScatter context on output
1521: */
1522: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1523: {
1524: PetscInt n, N;
1525: VecScatter vscat;
1526: Vec Y;
1528: PetscFunctionBegin;
1529: PetscCall(ISGetLocalSize(is, &n));
1530: PetscCall(ISGetSize(is, &N));
1531: PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1532: PetscCall(VecSetSizes(Y, n, N));
1533: PetscCall(VecSetBlockSize(Y, bs));
1534: PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1535: PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1536: PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1537: PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1538: PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1539: PetscCall(VecScatterDestroy(&vscat));
1540: *Z = Y;
1541: PetscFunctionReturn(PETSC_SUCCESS);
1542: }
1544: /*@
1545: VecGetSubVector - Gets a vector representing part of another vector
1547: Collective
1549: Input Parameters:
1550: + X - vector from which to extract a subvector
1551: - is - index set representing portion of `X` to extract
1553: Output Parameter:
1554: . Y - subvector corresponding to `is`
1556: Level: advanced
1558: Notes:
1559: The subvector `Y` should be returned with `VecRestoreSubVector()`.
1560: `X` and `is` must be defined on the same communicator
1562: Changes to the subvector will be reflected in the `X` vector on the call to `VecRestoreSubVector()`.
1564: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1565: modifying the subvector. Other non-overlapping subvectors can still be obtained from `X` using this function.
1567: 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`.
1569: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1570: @*/
1571: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1572: {
1573: Vec Z;
1575: PetscFunctionBegin;
1578: PetscCheckSameComm(X, 1, is, 2);
1579: PetscAssertPointer(Y, 3);
1580: if (X->ops->getsubvector) {
1581: PetscUseTypeMethod(X, getsubvector, is, &Z);
1582: } else { /* Default implementation currently does no caching */
1583: PetscBool contig;
1584: PetscInt n, N, start, bs;
1586: PetscCall(ISGetLocalSize(is, &n));
1587: PetscCall(ISGetSize(is, &N));
1588: PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1589: if (contig) { /* We can do a no-copy implementation */
1590: const PetscScalar *x;
1591: PetscInt state = 0;
1592: PetscBool isstd, iscuda, iship;
1594: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1595: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1596: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1597: if (iscuda) {
1598: #if defined(PETSC_HAVE_CUDA)
1599: const PetscScalar *x_d;
1600: PetscMPIInt size;
1601: PetscOffloadMask flg;
1603: PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1604: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1605: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1606: if (x) x += start;
1607: if (x_d) x_d += start;
1608: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1609: if (size == 1) {
1610: PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1611: } else {
1612: PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1613: }
1614: Z->offloadmask = flg;
1615: #endif
1616: } else if (iship) {
1617: #if defined(PETSC_HAVE_HIP)
1618: const PetscScalar *x_d;
1619: PetscMPIInt size;
1620: PetscOffloadMask flg;
1622: PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1623: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1624: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1625: if (x) x += start;
1626: if (x_d) x_d += start;
1627: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1628: if (size == 1) {
1629: PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1630: } else {
1631: PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1632: }
1633: Z->offloadmask = flg;
1634: #endif
1635: } else if (isstd) {
1636: PetscMPIInt size;
1638: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1639: PetscCall(VecGetArrayRead(X, &x));
1640: if (x) x += start;
1641: if (size == 1) {
1642: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1643: } else {
1644: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1645: }
1646: PetscCall(VecRestoreArrayRead(X, &x));
1647: } else { /* default implementation: use place array */
1648: PetscCall(VecGetArrayRead(X, &x));
1649: PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1650: PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1651: PetscCall(VecSetSizes(Z, n, N));
1652: PetscCall(VecSetBlockSize(Z, bs));
1653: PetscCall(VecPlaceArray(Z, PetscSafePointerPlusOffset(x, start)));
1654: PetscCall(VecRestoreArrayRead(X, &x));
1655: }
1657: /* this is relevant only in debug mode */
1658: PetscCall(VecLockGet(X, &state));
1659: if (state) PetscCall(VecLockReadPush(Z));
1660: Z->ops->placearray = NULL;
1661: Z->ops->replacearray = NULL;
1662: } else { /* Have to create a scatter and do a copy */
1663: PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1664: }
1665: }
1666: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1667: if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1668: PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1669: *Y = Z;
1670: PetscFunctionReturn(PETSC_SUCCESS);
1671: }
1673: /*@
1674: VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`
1676: Collective
1678: Input Parameters:
1679: + X - vector from which subvector was obtained
1680: . is - index set representing the subset of `X`
1681: - Y - subvector being restored
1683: Level: advanced
1685: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1686: @*/
1687: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1688: {
1689: PETSC_UNUSED PetscObjectState dummystate = 0;
1690: PetscBool unchanged;
1692: PetscFunctionBegin;
1695: PetscCheckSameComm(X, 1, is, 2);
1696: PetscAssertPointer(Y, 3);
1699: if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1700: else {
1701: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1702: if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1703: VecScatter scatter;
1704: PetscInt state;
1706: PetscCall(VecLockGet(X, &state));
1707: PetscCheck(state == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec X is locked for read-only or read/write access");
1709: PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1710: if (scatter) {
1711: PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1712: PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1713: } else {
1714: PetscBool iscuda, iship;
1715: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1716: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1718: if (iscuda) {
1719: #if defined(PETSC_HAVE_CUDA)
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(VecCUDAResetArray(*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(VecCUDAResetArray(*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 if (iship) {
1745: #if defined(PETSC_HAVE_HIP)
1746: PetscOffloadMask ymask = (*Y)->offloadmask;
1748: /* The offloadmask of X dictates where to move memory
1749: If X GPU data is valid, then move Y data on GPU if needed
1750: Otherwise, move back to the CPU */
1751: switch (X->offloadmask) {
1752: case PETSC_OFFLOAD_BOTH:
1753: if (ymask == PETSC_OFFLOAD_CPU) {
1754: PetscCall(VecHIPResetArray(*Y));
1755: } else if (ymask == PETSC_OFFLOAD_GPU) {
1756: X->offloadmask = PETSC_OFFLOAD_GPU;
1757: }
1758: break;
1759: case PETSC_OFFLOAD_GPU:
1760: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1761: break;
1762: case PETSC_OFFLOAD_CPU:
1763: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1764: break;
1765: case PETSC_OFFLOAD_UNALLOCATED:
1766: case PETSC_OFFLOAD_KOKKOS:
1767: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1768: }
1769: #endif
1770: } else {
1771: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1772: PetscCall(VecResetArray(*Y));
1773: }
1774: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1775: }
1776: }
1777: }
1778: PetscCall(VecDestroy(Y));
1779: PetscFunctionReturn(PETSC_SUCCESS);
1780: }
1782: /*@
1783: VecCreateLocalVector - Creates a vector object suitable for use with `VecGetLocalVector()` and friends. You must call `VecDestroy()` when the
1784: vector is no longer needed.
1786: Not Collective.
1788: Input Parameter:
1789: . v - The vector for which the local vector is desired.
1791: Output Parameter:
1792: . w - Upon exit this contains the local vector.
1794: Level: beginner
1796: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1797: @*/
1798: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1799: {
1800: VecType roottype;
1801: PetscInt n;
1803: PetscFunctionBegin;
1805: PetscAssertPointer(w, 2);
1806: if (v->ops->createlocalvector) {
1807: PetscUseTypeMethod(v, createlocalvector, w);
1808: PetscFunctionReturn(PETSC_SUCCESS);
1809: }
1810: PetscCall(VecGetRootType_Private(v, &roottype));
1811: PetscCall(VecCreate(PETSC_COMM_SELF, w));
1812: PetscCall(VecGetLocalSize(v, &n));
1813: PetscCall(VecSetSizes(*w, n, n));
1814: PetscCall(VecGetBlockSize(v, &n));
1815: PetscCall(VecSetBlockSize(*w, n));
1816: PetscCall(VecSetType(*w, roottype));
1817: PetscFunctionReturn(PETSC_SUCCESS);
1818: }
1820: /*@
1821: VecGetLocalVectorRead - Maps the local portion of a vector into a
1822: vector.
1824: Not Collective.
1826: Input Parameter:
1827: . v - The vector for which the local vector is desired.
1829: Output Parameter:
1830: . w - Upon exit this contains the local vector.
1832: Level: beginner
1834: Notes:
1835: You must call `VecRestoreLocalVectorRead()` when the local
1836: vector is no longer needed.
1838: This function is similar to `VecGetArrayRead()` which maps the local
1839: portion into a raw pointer. `VecGetLocalVectorRead()` is usually
1840: almost as efficient as `VecGetArrayRead()` but in certain circumstances
1841: `VecGetLocalVectorRead()` can be much more efficient than
1842: `VecGetArrayRead()`. This is because the construction of a contiguous
1843: array representing the vector data required by `VecGetArrayRead()` can
1844: be an expensive operation for certain vector types. For example, for
1845: GPU vectors `VecGetArrayRead()` requires that the data between device
1846: and host is synchronized.
1848: Unlike `VecGetLocalVector()`, this routine is not collective and
1849: preserves cached information.
1851: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1852: @*/
1853: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1854: {
1855: PetscFunctionBegin;
1858: VecCheckSameLocalSize(v, 1, w, 2);
1859: if (v->ops->getlocalvectorread) {
1860: PetscUseTypeMethod(v, getlocalvectorread, w);
1861: } else {
1862: PetscScalar *a;
1864: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1865: PetscCall(VecPlaceArray(w, a));
1866: }
1867: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1868: PetscCall(VecLockReadPush(v));
1869: PetscCall(VecLockReadPush(w));
1870: PetscFunctionReturn(PETSC_SUCCESS);
1871: }
1873: /*@
1874: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1875: previously mapped into a vector using `VecGetLocalVectorRead()`.
1877: Not Collective.
1879: Input Parameters:
1880: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
1881: - w - The vector into which the local portion of `v` was mapped.
1883: Level: beginner
1885: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1886: @*/
1887: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1888: {
1889: PetscFunctionBegin;
1892: if (v->ops->restorelocalvectorread) {
1893: PetscUseTypeMethod(v, restorelocalvectorread, w);
1894: } else {
1895: const PetscScalar *a;
1897: PetscCall(VecGetArrayRead(w, &a));
1898: PetscCall(VecRestoreArrayRead(v, &a));
1899: PetscCall(VecResetArray(w));
1900: }
1901: PetscCall(VecLockReadPop(v));
1902: PetscCall(VecLockReadPop(w));
1903: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1904: PetscFunctionReturn(PETSC_SUCCESS);
1905: }
1907: /*@
1908: VecGetLocalVector - Maps the local portion of a vector into a
1909: vector.
1911: Collective
1913: Input Parameter:
1914: . v - The vector for which the local vector is desired.
1916: Output Parameter:
1917: . w - Upon exit this contains the local vector.
1919: Level: beginner
1921: Notes:
1922: You must call `VecRestoreLocalVector()` when the local
1923: vector is no longer needed.
1925: This function is similar to `VecGetArray()` which maps the local
1926: portion into a raw pointer. `VecGetLocalVector()` is usually about as
1927: efficient as `VecGetArray()` but in certain circumstances
1928: `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1929: This is because the construction of a contiguous array representing
1930: the vector data required by `VecGetArray()` can be an expensive
1931: operation for certain vector types. For example, for GPU vectors
1932: `VecGetArray()` requires that the data between device and host is
1933: synchronized.
1935: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1936: @*/
1937: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1938: {
1939: PetscFunctionBegin;
1942: VecCheckSameLocalSize(v, 1, w, 2);
1943: if (v->ops->getlocalvector) {
1944: PetscUseTypeMethod(v, getlocalvector, w);
1945: } else {
1946: PetscScalar *a;
1948: PetscCall(VecGetArray(v, &a));
1949: PetscCall(VecPlaceArray(w, a));
1950: }
1951: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1952: PetscFunctionReturn(PETSC_SUCCESS);
1953: }
1955: /*@
1956: VecRestoreLocalVector - Unmaps the local portion of a vector
1957: previously mapped into a vector using `VecGetLocalVector()`.
1959: Logically Collective.
1961: Input Parameters:
1962: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
1963: - w - The vector into which the local portion of `v` was mapped.
1965: Level: beginner
1967: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1968: @*/
1969: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1970: {
1971: PetscFunctionBegin;
1974: if (v->ops->restorelocalvector) {
1975: PetscUseTypeMethod(v, restorelocalvector, w);
1976: } else {
1977: PetscScalar *a;
1978: PetscCall(VecGetArray(w, &a));
1979: PetscCall(VecRestoreArray(v, &a));
1980: PetscCall(VecResetArray(w));
1981: }
1982: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1983: PetscCall(PetscObjectStateIncrease((PetscObject)v));
1984: PetscFunctionReturn(PETSC_SUCCESS);
1985: }
1987: /*@C
1988: VecGetArray - Returns a pointer to a contiguous array that contains this
1989: MPI processes's portion of the vector data
1991: Logically Collective
1993: Input Parameter:
1994: . x - the vector
1996: Output Parameter:
1997: . a - location to put pointer to the array
1999: Level: beginner
2001: Notes:
2002: For the standard PETSc vectors, `VecGetArray()` returns a pointer to the local data array and
2003: does not use any copies. If the underlying vector data is not stored in a contiguous array
2004: this routine will copy the data to a contiguous array and return a pointer to that. You MUST
2005: call `VecRestoreArray()` when you no longer need access to the array.
2007: Fortran Notes:
2008: `VecGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayF90()`
2010: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
2011: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2012: @*/
2013: PetscErrorCode VecGetArray(Vec x, PetscScalar **a)
2014: {
2015: PetscFunctionBegin;
2017: PetscCall(VecSetErrorIfLocked(x, 1));
2018: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
2019: PetscUseTypeMethod(x, getarray, a);
2020: } else if (x->petscnative) { /* VECSTANDARD */
2021: *a = *((PetscScalar **)x->data);
2022: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
2023: PetscFunctionReturn(PETSC_SUCCESS);
2024: }
2026: /*@C
2027: VecRestoreArray - Restores a vector after `VecGetArray()` has been called and the array is no longer needed
2029: Logically Collective
2031: Input Parameters:
2032: + x - the vector
2033: - a - location of pointer to array obtained from `VecGetArray()`
2035: Level: beginner
2037: Fortran Notes:
2038: `VecRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayF90()`
2040: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2041: `VecGetArrayPair()`, `VecRestoreArrayPair()`
2042: @*/
2043: PetscErrorCode VecRestoreArray(Vec x, PetscScalar **a)
2044: {
2045: PetscFunctionBegin;
2047: if (a) PetscAssertPointer(a, 2);
2048: if (x->ops->restorearray) {
2049: PetscUseTypeMethod(x, restorearray, a);
2050: } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
2051: if (a) *a = NULL;
2052: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2053: PetscFunctionReturn(PETSC_SUCCESS);
2054: }
2055: /*@C
2056: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
2058: Not Collective
2060: Input Parameter:
2061: . x - the vector
2063: Output Parameter:
2064: . a - the array
2066: Level: beginner
2068: Notes:
2069: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2071: Unlike `VecGetArray()`, preserves cached information like vector norms.
2073: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
2074: implementations may require a copy, but such implementations should cache the contiguous representation so that
2075: only one copy is performed when this routine is called multiple times in sequence.
2077: Fortran Notes:
2078: `VecGetArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayReadF90()`
2080: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2081: @*/
2082: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar **a)
2083: {
2084: PetscFunctionBegin;
2086: PetscAssertPointer(a, 2);
2087: if (x->ops->getarrayread) {
2088: PetscUseTypeMethod(x, getarrayread, a);
2089: } else if (x->ops->getarray) {
2090: PetscObjectState state;
2092: /* VECNEST, VECCUDA, VECKOKKOS etc */
2093: // x->ops->getarray may bump the object state, but since we know this is a read-only get
2094: // we can just undo that
2095: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2096: PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
2097: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2098: } else if (x->petscnative) {
2099: /* VECSTANDARD */
2100: *a = *((PetscScalar **)x->data);
2101: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
2102: PetscFunctionReturn(PETSC_SUCCESS);
2103: }
2105: /*@C
2106: VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`
2108: Not Collective
2110: Input Parameters:
2111: + x - the vector
2112: - a - the array
2114: Level: beginner
2116: Fortran Notes:
2117: `VecRestoreArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayReadF90()`
2119: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2120: @*/
2121: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar **a)
2122: {
2123: PetscFunctionBegin;
2125: if (a) PetscAssertPointer(a, 2);
2126: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
2127: /* nothing */
2128: } else if (x->ops->restorearrayread) { /* VECNEST */
2129: PetscUseTypeMethod(x, restorearrayread, a);
2130: } else { /* No one? */
2131: PetscObjectState state;
2133: // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2134: // we can just undo that
2135: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2136: PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2137: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2138: }
2139: if (a) *a = NULL;
2140: PetscFunctionReturn(PETSC_SUCCESS);
2141: }
2143: /*@C
2144: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
2145: MPI processes's portion of the vector data.
2147: Logically Collective
2149: Input Parameter:
2150: . x - the vector
2152: Output Parameter:
2153: . a - location to put pointer to the array
2155: Level: intermediate
2157: Note:
2158: The values in this array are NOT valid, the caller of this routine is responsible for putting
2159: values into the array; any values it does not set will be invalid.
2161: The array must be returned using a matching call to `VecRestoreArrayWrite()`.
2163: For vectors associated with GPUs, the host and device vectors are not synchronized before
2164: giving access. If you need correct values in the array use `VecGetArray()`
2166: Fortran Notes:
2167: `VecGetArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayF90()`
2169: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
2170: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
2171: @*/
2172: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar **a)
2173: {
2174: PetscFunctionBegin;
2176: PetscAssertPointer(a, 2);
2177: PetscCall(VecSetErrorIfLocked(x, 1));
2178: if (x->ops->getarraywrite) {
2179: PetscUseTypeMethod(x, getarraywrite, a);
2180: } else {
2181: PetscCall(VecGetArray(x, a));
2182: }
2183: PetscFunctionReturn(PETSC_SUCCESS);
2184: }
2186: /*@C
2187: VecRestoreArrayWrite - Restores a vector after `VecGetArrayWrite()` has been called.
2189: Logically Collective
2191: Input Parameters:
2192: + x - the vector
2193: - a - location of pointer to array obtained from `VecGetArray()`
2195: Level: beginner
2197: Fortran Notes:
2198: `VecRestoreArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayF90()`
2200: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2201: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2202: @*/
2203: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar **a)
2204: {
2205: PetscFunctionBegin;
2207: if (a) PetscAssertPointer(a, 2);
2208: if (x->ops->restorearraywrite) {
2209: PetscUseTypeMethod(x, restorearraywrite, a);
2210: } else if (x->ops->restorearray) {
2211: PetscUseTypeMethod(x, restorearray, a);
2212: }
2213: if (a) *a = NULL;
2214: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2215: PetscFunctionReturn(PETSC_SUCCESS);
2216: }
2218: /*@C
2219: VecGetArrays - Returns a pointer to the arrays in a set of vectors
2220: that were created by a call to `VecDuplicateVecs()`.
2222: Logically Collective; No Fortran Support
2224: Input Parameters:
2225: + x - the vectors
2226: - n - the number of vectors
2228: Output Parameter:
2229: . a - location to put pointer to the array
2231: Level: intermediate
2233: Note:
2234: You MUST call `VecRestoreArrays()` when you no longer need access to the arrays.
2236: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2237: @*/
2238: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2239: {
2240: PetscInt i;
2241: PetscScalar **q;
2243: PetscFunctionBegin;
2244: PetscAssertPointer(x, 1);
2246: PetscAssertPointer(a, 3);
2247: PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2248: PetscCall(PetscMalloc1(n, &q));
2249: for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2250: *a = q;
2251: PetscFunctionReturn(PETSC_SUCCESS);
2252: }
2254: /*@C
2255: VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2256: has been called.
2258: Logically Collective; No Fortran Support
2260: Input Parameters:
2261: + x - the vector
2262: . n - the number of vectors
2263: - a - location of pointer to arrays obtained from `VecGetArrays()`
2265: Notes:
2266: For regular PETSc vectors this routine does not involve any copies. For
2267: any special vectors that do not store local vector data in a contiguous
2268: array, this routine will copy the data back into the underlying
2269: vector data structure from the arrays obtained with `VecGetArrays()`.
2271: Level: intermediate
2273: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2274: @*/
2275: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2276: {
2277: PetscInt i;
2278: PetscScalar **q = *a;
2280: PetscFunctionBegin;
2281: PetscAssertPointer(x, 1);
2283: PetscAssertPointer(a, 3);
2285: for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2286: PetscCall(PetscFree(q));
2287: PetscFunctionReturn(PETSC_SUCCESS);
2288: }
2290: /*@C
2291: VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g.,
2292: `VECCUDA`), the returned pointer will be a device pointer to the device memory that contains
2293: this MPI processes's portion of the vector data.
2295: Logically Collective; No Fortran Support
2297: Input Parameter:
2298: . x - the vector
2300: Output Parameters:
2301: + a - location to put pointer to the array
2302: - mtype - memory type of the array
2304: Level: beginner
2306: Note:
2307: Device data is guaranteed to have the latest value. Otherwise, when this is a host vector
2308: (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host
2309: pointer.
2311: For `VECKOKKOS`, if Kokkos is configured without device (e.g., use serial or openmp), per
2312: this function, the vector works like `VECSEQ`/`VECMPI`; otherwise, it works like `VECCUDA` or
2313: `VECHIP` etc.
2315: Use `VecRestoreArrayAndMemType()` when the array access is no longer needed.
2317: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`,
2318: `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2319: @*/
2320: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2321: {
2322: PetscFunctionBegin;
2325: PetscAssertPointer(a, 2);
2326: if (mtype) PetscAssertPointer(mtype, 3);
2327: PetscCall(VecSetErrorIfLocked(x, 1));
2328: if (x->ops->getarrayandmemtype) {
2329: /* VECCUDA, VECKOKKOS etc */
2330: PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2331: } else {
2332: /* VECSTANDARD, VECNEST, VECVIENNACL */
2333: PetscCall(VecGetArray(x, a));
2334: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2335: }
2336: PetscFunctionReturn(PETSC_SUCCESS);
2337: }
2339: /*@C
2340: VecRestoreArrayAndMemType - Restores a vector after `VecGetArrayAndMemType()` has been called.
2342: Logically Collective; No Fortran Support
2344: Input Parameters:
2345: + x - the vector
2346: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`
2348: Level: beginner
2350: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`,
2351: `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2352: @*/
2353: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar **a)
2354: {
2355: PetscFunctionBegin;
2358: if (a) PetscAssertPointer(a, 2);
2359: if (x->ops->restorearrayandmemtype) {
2360: /* VECCUDA, VECKOKKOS etc */
2361: PetscUseTypeMethod(x, restorearrayandmemtype, a);
2362: } else {
2363: /* VECNEST, VECVIENNACL */
2364: PetscCall(VecRestoreArray(x, a));
2365: } /* VECSTANDARD does nothing */
2366: if (a) *a = NULL;
2367: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2368: PetscFunctionReturn(PETSC_SUCCESS);
2369: }
2371: /*@C
2372: VecGetArrayReadAndMemType - Like `VecGetArrayRead()`, but if the input vector is a device vector, it will return a read-only device pointer.
2373: The returned pointer is guaranteed to point to up-to-date data. For host vectors, it functions as `VecGetArrayRead()`.
2375: Not Collective; No Fortran Support
2377: Input Parameter:
2378: . x - the vector
2380: Output Parameters:
2381: + a - the array
2382: - mtype - memory type of the array
2384: Level: beginner
2386: Notes:
2387: The array must be returned using a matching call to `VecRestoreArrayReadAndMemType()`.
2389: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2390: @*/
2391: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar **a, PetscMemType *mtype)
2392: {
2393: PetscFunctionBegin;
2396: PetscAssertPointer(a, 2);
2397: if (mtype) PetscAssertPointer(mtype, 3);
2398: if (x->ops->getarrayreadandmemtype) {
2399: /* VECCUDA/VECHIP though they are also petscnative */
2400: PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2401: } else if (x->ops->getarrayandmemtype) {
2402: /* VECKOKKOS */
2403: PetscObjectState state;
2405: // see VecGetArrayRead() for why
2406: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2407: PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2408: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2409: } else {
2410: PetscCall(VecGetArrayRead(x, a));
2411: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2412: }
2413: PetscFunctionReturn(PETSC_SUCCESS);
2414: }
2416: /*@C
2417: VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`
2419: Not Collective; No Fortran Support
2421: Input Parameters:
2422: + x - the vector
2423: - a - the array
2425: Level: beginner
2427: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2428: @*/
2429: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar **a)
2430: {
2431: PetscFunctionBegin;
2434: if (a) PetscAssertPointer(a, 2);
2435: if (x->ops->restorearrayreadandmemtype) {
2436: /* VECCUDA/VECHIP */
2437: PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2438: } else if (!x->petscnative) {
2439: /* VECNEST */
2440: PetscCall(VecRestoreArrayRead(x, a));
2441: }
2442: if (a) *a = NULL;
2443: PetscFunctionReturn(PETSC_SUCCESS);
2444: }
2446: /*@C
2447: VecGetArrayWriteAndMemType - Like `VecGetArrayWrite()`, but if this is a device vector it will always return
2448: a device pointer to the device memory that contains this processor's portion of the vector data.
2450: Logically Collective; No Fortran Support
2452: Input Parameter:
2453: . x - the vector
2455: Output Parameters:
2456: + a - the array
2457: - mtype - memory type of the array
2459: Level: beginner
2461: Note:
2462: The array must be returned using a matching call to `VecRestoreArrayWriteAndMemType()`, where it will label the device memory as most recent.
2464: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2465: @*/
2466: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2467: {
2468: PetscFunctionBegin;
2471: PetscCall(VecSetErrorIfLocked(x, 1));
2472: PetscAssertPointer(a, 2);
2473: if (mtype) PetscAssertPointer(mtype, 3);
2474: if (x->ops->getarraywriteandmemtype) {
2475: /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2476: PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2477: } else if (x->ops->getarrayandmemtype) {
2478: PetscCall(VecGetArrayAndMemType(x, a, mtype));
2479: } else {
2480: /* VECNEST, VECVIENNACL */
2481: PetscCall(VecGetArrayWrite(x, a));
2482: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2483: }
2484: PetscFunctionReturn(PETSC_SUCCESS);
2485: }
2487: /*@C
2488: VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`
2490: Logically Collective; No Fortran Support
2492: Input Parameters:
2493: + x - the vector
2494: - a - the array
2496: Level: beginner
2498: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2499: @*/
2500: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar **a)
2501: {
2502: PetscFunctionBegin;
2505: PetscCall(VecSetErrorIfLocked(x, 1));
2506: if (a) PetscAssertPointer(a, 2);
2507: if (x->ops->restorearraywriteandmemtype) {
2508: /* VECCUDA/VECHIP */
2509: PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2510: PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2511: } else if (x->ops->restorearrayandmemtype) {
2512: PetscCall(VecRestoreArrayAndMemType(x, a));
2513: } else {
2514: PetscCall(VecRestoreArray(x, a));
2515: }
2516: if (a) *a = NULL;
2517: PetscFunctionReturn(PETSC_SUCCESS);
2518: }
2520: /*@
2521: VecPlaceArray - Allows one to replace the array in a vector with an
2522: array provided by the user. This is useful to avoid copying an array
2523: into a vector.
2525: Logically Collective; No Fortran Support
2527: Input Parameters:
2528: + vec - the vector
2529: - array - the array
2531: Level: developer
2533: Notes:
2534: Adding `const` to `array` was an oversight, as subsequent operations on `vec` would
2535: likely modify the data in `array`. However, we have kept it to avoid breaking APIs.
2537: Use `VecReplaceArray()` instead to permanently replace the array
2539: You can return to the original array with a call to `VecResetArray()`. `vec` does not take
2540: ownership of `array` in any way.
2542: The user must free `array` themselves but be careful not to
2543: do so before the vector has either been destroyed, had its original array restored with
2544: `VecResetArray()` or permanently replaced with `VecReplaceArray()`.
2546: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2547: @*/
2548: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2549: {
2550: PetscFunctionBegin;
2553: if (array) PetscAssertPointer(array, 2);
2554: PetscUseTypeMethod(vec, placearray, array);
2555: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2556: PetscFunctionReturn(PETSC_SUCCESS);
2557: }
2559: /*@C
2560: VecReplaceArray - Allows one to replace the array in a vector with an
2561: array provided by the user. This is useful to avoid copying an array
2562: into a vector.
2564: Logically Collective; No Fortran Support
2566: Input Parameters:
2567: + vec - the vector
2568: - array - the array
2570: Level: developer
2572: Notes:
2573: Adding `const` to `array` was an oversight, as subsequent operations on `vec` would
2574: likely modify the data in `array`. However, we have kept it to avoid breaking APIs.
2576: This permanently replaces the array and frees the memory associated
2577: with the old array. Use `VecPlaceArray()` to temporarily replace the array.
2579: The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
2580: freed by the user. It will be freed when the vector is destroyed.
2582: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2583: @*/
2584: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2585: {
2586: PetscFunctionBegin;
2589: PetscUseTypeMethod(vec, replacearray, array);
2590: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2591: PetscFunctionReturn(PETSC_SUCCESS);
2592: }
2594: /*MC
2595: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
2596: and makes them accessible via a Fortran pointer.
2598: Synopsis:
2599: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
2601: Collective
2603: Input Parameters:
2604: + x - a vector to mimic
2605: - n - the number of vectors to obtain
2607: Output Parameters:
2608: + y - Fortran pointer to the array of vectors
2609: - ierr - error code
2611: Example of Usage:
2612: .vb
2613: #include <petsc/finclude/petscvec.h>
2614: use petscvec
2616: Vec x
2617: Vec, pointer :: y(:)
2618: ....
2619: call VecDuplicateVecsF90(x,2,y,ierr)
2620: call VecSet(y(2),alpha,ierr)
2621: call VecSet(y(2),alpha,ierr)
2622: ....
2623: call VecDestroyVecsF90(2,y,ierr)
2624: .ve
2626: Level: beginner
2628: Note:
2629: Use `VecDestroyVecsF90()` to free the space.
2631: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecsF90()`, `VecDuplicateVecs()`
2632: M*/
2634: /*MC
2635: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2636: `VecGetArrayF90()`.
2638: Synopsis:
2639: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2641: Logically Collective
2643: Input Parameters:
2644: + x - vector
2645: - xx_v - the Fortran pointer to the array
2647: Output Parameter:
2648: . ierr - error code
2650: Example of Usage:
2651: .vb
2652: #include <petsc/finclude/petscvec.h>
2653: use petscvec
2655: PetscScalar, pointer :: xx_v(:)
2656: ....
2657: call VecGetArrayF90(x,xx_v,ierr)
2658: xx_v(3) = a
2659: call VecRestoreArrayF90(x,xx_v,ierr)
2660: .ve
2662: Level: beginner
2664: .seealso: [](ch_vectors), `Vec`, `VecGetArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrayReadF90()`
2665: M*/
2667: /*MC
2668: VecDestroyVecsF90 - Frees a block of vectors obtained with `VecDuplicateVecsF90()`.
2670: Synopsis:
2671: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
2673: Collective
2675: Input Parameters:
2676: + n - the number of vectors previously obtained
2677: - x - pointer to array of vector pointers
2679: Output Parameter:
2680: . ierr - error code
2682: Level: beginner
2684: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecs()`, `VecDuplicateVecsF90()`
2685: M*/
2687: /*MC
2688: VecGetArrayF90 - Accesses a vector array from Fortran. For default PETSc
2689: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2690: this routine is implementation dependent. You MUST call `VecRestoreArrayF90()`
2691: when you no longer need access to the array.
2693: Synopsis:
2694: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2696: Logically Collective
2698: Input Parameter:
2699: . x - vector
2701: Output Parameters:
2702: + xx_v - the Fortran pointer to the array
2703: - ierr - error code
2705: Example of Usage:
2706: .vb
2707: #include <petsc/finclude/petscvec.h>
2708: use petscvec
2710: PetscScalar, pointer :: xx_v(:)
2711: ....
2712: call VecGetArrayF90(x,xx_v,ierr)
2713: xx_v(3) = a
2714: call VecRestoreArrayF90(x,xx_v,ierr)
2715: .ve
2717: Level: beginner
2719: Note:
2720: If you ONLY intend to read entries from the array and not change any entries you should use `VecGetArrayReadF90()`.
2722: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayReadF90()`
2723: M*/
2725: /*MC
2726: VecGetArrayReadF90 - Accesses a read only array from Fortran. For default PETSc
2727: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2728: this routine is implementation dependent. You MUST call `VecRestoreArrayReadF90()`
2729: when you no longer need access to the array.
2731: Synopsis:
2732: VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2734: Logically Collective
2736: Input Parameter:
2737: . x - vector
2739: Output Parameters:
2740: + xx_v - the Fortran pointer to the array
2741: - ierr - error code
2743: Example of Usage:
2744: .vb
2745: #include <petsc/finclude/petscvec.h>
2746: use petscvec
2748: PetscScalar, pointer :: xx_v(:)
2749: ....
2750: call VecGetArrayReadF90(x,xx_v,ierr)
2751: a = xx_v(3)
2752: call VecRestoreArrayReadF90(x,xx_v,ierr)
2753: .ve
2755: Level: beginner
2757: Note:
2758: If you intend to write entries into the array you must use `VecGetArrayF90()`.
2760: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecGetArrayF90()`
2761: M*/
2763: /*MC
2764: VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2765: `VecGetArrayReadF90()`.
2767: Synopsis:
2768: VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2770: Logically Collective
2772: Input Parameters:
2773: + x - vector
2774: - xx_v - the Fortran pointer to the array
2776: Output Parameter:
2777: . ierr - error code
2779: Example of Usage:
2780: .vb
2781: #include <petsc/finclude/petscvec.h>
2782: use petscvec
2784: PetscScalar, pointer :: xx_v(:)
2785: ....
2786: call VecGetArrayReadF90(x,xx_v,ierr)
2787: a = xx_v(3)
2788: call VecRestoreArrayReadF90(x,xx_v,ierr)
2789: .ve
2791: Level: beginner
2793: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecRestoreArrayF90()`
2794: M*/
2796: /*@C
2797: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2798: processor's portion of the vector data. You MUST call `VecRestoreArray2d()`
2799: when you no longer need access to the array.
2801: Logically Collective
2803: Input Parameters:
2804: + x - the vector
2805: . m - first dimension of two dimensional array
2806: . n - second dimension of two dimensional array
2807: . mstart - first index you will use in first coordinate direction (often 0)
2808: - nstart - first index in the second coordinate direction (often 0)
2810: Output Parameter:
2811: . a - location to put pointer to the array
2813: Level: developer
2815: Notes:
2816: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2817: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2818: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2819: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2821: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2823: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2824: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2825: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2826: @*/
2827: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2828: {
2829: PetscInt i, N;
2830: PetscScalar *aa;
2832: PetscFunctionBegin;
2834: PetscAssertPointer(a, 6);
2836: PetscCall(VecGetLocalSize(x, &N));
2837: PetscCheck(m * n == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local array size %" PetscInt_FMT " does not match 2d array dimensions %" PetscInt_FMT " by %" PetscInt_FMT, N, m, n);
2838: PetscCall(VecGetArray(x, &aa));
2840: PetscCall(PetscMalloc1(m, a));
2841: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2842: *a -= mstart;
2843: PetscFunctionReturn(PETSC_SUCCESS);
2844: }
2846: /*@C
2847: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2848: processor's portion of the vector data. You MUST call `VecRestoreArray2dWrite()`
2849: when you no longer need access to the array.
2851: Logically Collective
2853: Input Parameters:
2854: + x - the vector
2855: . m - first dimension of two dimensional array
2856: . n - second dimension of two dimensional array
2857: . mstart - first index you will use in first coordinate direction (often 0)
2858: - nstart - first index in the second coordinate direction (often 0)
2860: Output Parameter:
2861: . a - location to put pointer to the array
2863: Level: developer
2865: Notes:
2866: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2867: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2868: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2869: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2871: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2873: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2874: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2875: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2876: @*/
2877: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2878: {
2879: PetscInt i, N;
2880: PetscScalar *aa;
2882: PetscFunctionBegin;
2884: PetscAssertPointer(a, 6);
2886: PetscCall(VecGetLocalSize(x, &N));
2887: 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);
2888: PetscCall(VecGetArrayWrite(x, &aa));
2890: PetscCall(PetscMalloc1(m, a));
2891: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2892: *a -= mstart;
2893: PetscFunctionReturn(PETSC_SUCCESS);
2894: }
2896: /*@C
2897: VecRestoreArray2d - Restores a vector after `VecGetArray2d()` has been called.
2899: Logically Collective
2901: Input Parameters:
2902: + x - the vector
2903: . m - first dimension of two dimensional array
2904: . n - second dimension of the two dimensional array
2905: . mstart - first index you will use in first coordinate direction (often 0)
2906: . nstart - first index in the second coordinate direction (often 0)
2907: - a - location of pointer to array obtained from `VecGetArray2d()`
2909: Level: developer
2911: Notes:
2912: For regular PETSc vectors this routine does not involve any copies. For
2913: any special vectors that do not store local vector data in a contiguous
2914: array, this routine will copy the data back into the underlying
2915: vector data structure from the array obtained with `VecGetArray()`.
2917: This routine actually zeros out the `a` pointer.
2919: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2920: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2921: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2922: @*/
2923: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2924: {
2925: void *dummy;
2927: PetscFunctionBegin;
2929: PetscAssertPointer(a, 6);
2931: dummy = (void *)(*a + mstart);
2932: PetscCall(PetscFree(dummy));
2933: PetscCall(VecRestoreArray(x, NULL));
2934: *a = NULL;
2935: PetscFunctionReturn(PETSC_SUCCESS);
2936: }
2938: /*@C
2939: VecRestoreArray2dWrite - Restores a vector after `VecGetArray2dWrite()` has been called.
2941: Logically Collective
2943: Input Parameters:
2944: + x - the vector
2945: . m - first dimension of two dimensional array
2946: . n - second dimension of the two dimensional array
2947: . mstart - first index you will use in first coordinate direction (often 0)
2948: . nstart - first index in the second coordinate direction (often 0)
2949: - a - location of pointer to array obtained from `VecGetArray2d()`
2951: Level: developer
2953: Notes:
2954: For regular PETSc vectors this routine does not involve any copies. For
2955: any special vectors that do not store local vector data in a contiguous
2956: array, this routine will copy the data back into the underlying
2957: vector data structure from the array obtained with `VecGetArray()`.
2959: This routine actually zeros out the `a` pointer.
2961: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2962: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2963: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2964: @*/
2965: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2966: {
2967: void *dummy;
2969: PetscFunctionBegin;
2971: PetscAssertPointer(a, 6);
2973: dummy = (void *)(*a + mstart);
2974: PetscCall(PetscFree(dummy));
2975: PetscCall(VecRestoreArrayWrite(x, NULL));
2976: PetscFunctionReturn(PETSC_SUCCESS);
2977: }
2979: /*@C
2980: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2981: processor's portion of the vector data. You MUST call `VecRestoreArray1d()`
2982: when you no longer need access to the array.
2984: Logically Collective
2986: Input Parameters:
2987: + x - the vector
2988: . m - first dimension of two dimensional array
2989: - mstart - first index you will use in first coordinate direction (often 0)
2991: Output Parameter:
2992: . a - location to put pointer to the array
2994: Level: developer
2996: Notes:
2997: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2998: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2999: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3001: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3003: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3004: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3005: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3006: @*/
3007: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3008: {
3009: PetscInt N;
3011: PetscFunctionBegin;
3013: PetscAssertPointer(a, 4);
3015: PetscCall(VecGetLocalSize(x, &N));
3016: PetscCheck(m == N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local array size %" PetscInt_FMT " does not match 1d array dimensions %" PetscInt_FMT, N, m);
3017: PetscCall(VecGetArray(x, a));
3018: *a -= mstart;
3019: PetscFunctionReturn(PETSC_SUCCESS);
3020: }
3022: /*@C
3023: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
3024: processor's portion of the vector data. You MUST call `VecRestoreArray1dWrite()`
3025: when you no longer need access to the array.
3027: Logically Collective
3029: Input Parameters:
3030: + x - the vector
3031: . m - first dimension of two dimensional array
3032: - mstart - first index you will use in first coordinate direction (often 0)
3034: Output Parameter:
3035: . a - location to put pointer to the array
3037: Level: developer
3039: Notes:
3040: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3041: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3042: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3044: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3046: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3047: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3048: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3049: @*/
3050: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3051: {
3052: PetscInt N;
3054: PetscFunctionBegin;
3056: PetscAssertPointer(a, 4);
3058: PetscCall(VecGetLocalSize(x, &N));
3059: 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);
3060: PetscCall(VecGetArrayWrite(x, a));
3061: *a -= mstart;
3062: PetscFunctionReturn(PETSC_SUCCESS);
3063: }
3065: /*@C
3066: VecRestoreArray1d - Restores a vector after `VecGetArray1d()` has been called.
3068: Logically Collective
3070: Input Parameters:
3071: + x - the vector
3072: . m - first dimension of two dimensional array
3073: . mstart - first index you will use in first coordinate direction (often 0)
3074: - a - location of pointer to array obtained from `VecGetArray1d()`
3076: Level: developer
3078: Notes:
3079: For regular PETSc vectors this routine does not involve any copies. For
3080: any special vectors that do not store local vector data in a contiguous
3081: array, this routine will copy the data back into the underlying
3082: vector data structure from the array obtained with `VecGetArray1d()`.
3084: This routine actually zeros out the `a` pointer.
3086: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3087: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3088: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3089: @*/
3090: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3091: {
3092: PetscFunctionBegin;
3095: PetscCall(VecRestoreArray(x, NULL));
3096: *a = NULL;
3097: PetscFunctionReturn(PETSC_SUCCESS);
3098: }
3100: /*@C
3101: VecRestoreArray1dWrite - Restores a vector after `VecGetArray1dWrite()` has been called.
3103: Logically Collective
3105: Input Parameters:
3106: + x - the vector
3107: . m - first dimension of two dimensional array
3108: . mstart - first index you will use in first coordinate direction (often 0)
3109: - a - location of pointer to array obtained from `VecGetArray1d()`
3111: Level: developer
3113: Notes:
3114: For regular PETSc vectors this routine does not involve any copies. For
3115: any special vectors that do not store local vector data in a contiguous
3116: array, this routine will copy the data back into the underlying
3117: vector data structure from the array obtained with `VecGetArray1d()`.
3119: This routine actually zeros out the `a` pointer.
3121: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3122: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3123: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3124: @*/
3125: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3126: {
3127: PetscFunctionBegin;
3130: PetscCall(VecRestoreArrayWrite(x, NULL));
3131: *a = NULL;
3132: PetscFunctionReturn(PETSC_SUCCESS);
3133: }
3135: /*@C
3136: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
3137: processor's portion of the vector data. You MUST call `VecRestoreArray3d()`
3138: when you no longer need access to the array.
3140: Logically Collective
3142: Input Parameters:
3143: + x - the vector
3144: . m - first dimension of three dimensional array
3145: . n - second dimension of three dimensional array
3146: . p - third dimension of three dimensional array
3147: . mstart - first index you will use in first coordinate direction (often 0)
3148: . nstart - first index in the second coordinate direction (often 0)
3149: - pstart - first index in the third coordinate direction (often 0)
3151: Output Parameter:
3152: . a - location to put pointer to the array
3154: Level: developer
3156: Notes:
3157: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3158: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3159: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3160: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3162: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3164: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3165: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
3166: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3167: @*/
3168: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3169: {
3170: PetscInt i, N, j;
3171: PetscScalar *aa, **b;
3173: PetscFunctionBegin;
3175: PetscAssertPointer(a, 8);
3177: PetscCall(VecGetLocalSize(x, &N));
3178: 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);
3179: PetscCall(VecGetArray(x, &aa));
3181: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3182: b = (PetscScalar **)((*a) + m);
3183: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3184: for (i = 0; i < m; i++)
3185: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset(aa, i * n * p + j * p - pstart);
3186: *a -= mstart;
3187: PetscFunctionReturn(PETSC_SUCCESS);
3188: }
3190: /*@C
3191: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
3192: processor's portion of the vector data. You MUST call `VecRestoreArray3dWrite()`
3193: when you no longer need access to the array.
3195: Logically Collective
3197: Input Parameters:
3198: + x - the vector
3199: . m - first dimension of three dimensional array
3200: . n - second dimension of three dimensional array
3201: . p - third dimension of three dimensional array
3202: . mstart - first index you will use in first coordinate direction (often 0)
3203: . nstart - first index in the second coordinate direction (often 0)
3204: - pstart - first index in the third coordinate direction (often 0)
3206: Output Parameter:
3207: . a - location to put pointer to the array
3209: Level: developer
3211: Notes:
3212: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3213: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3214: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3215: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3217: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3219: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3220: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3221: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3222: @*/
3223: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3224: {
3225: PetscInt i, N, j;
3226: PetscScalar *aa, **b;
3228: PetscFunctionBegin;
3230: PetscAssertPointer(a, 8);
3232: PetscCall(VecGetLocalSize(x, &N));
3233: 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);
3234: PetscCall(VecGetArrayWrite(x, &aa));
3236: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3237: b = (PetscScalar **)((*a) + m);
3238: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3239: for (i = 0; i < m; i++)
3240: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3242: *a -= mstart;
3243: PetscFunctionReturn(PETSC_SUCCESS);
3244: }
3246: /*@C
3247: VecRestoreArray3d - Restores a vector after `VecGetArray3d()` has been called.
3249: Logically Collective
3251: Input Parameters:
3252: + x - the vector
3253: . m - first dimension of three dimensional array
3254: . n - second dimension of the three dimensional array
3255: . p - third dimension of the three dimensional array
3256: . mstart - first index you will use in first coordinate direction (often 0)
3257: . nstart - first index in the second coordinate direction (often 0)
3258: . pstart - first index in the third coordinate direction (often 0)
3259: - a - location of pointer to array obtained from VecGetArray3d()
3261: Level: developer
3263: Notes:
3264: For regular PETSc vectors this routine does not involve any copies. For
3265: any special vectors that do not store local vector data in a contiguous
3266: array, this routine will copy the data back into the underlying
3267: vector data structure from the array obtained with `VecGetArray()`.
3269: This routine actually zeros out the `a` pointer.
3271: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3272: `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3273: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3274: @*/
3275: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3276: {
3277: void *dummy;
3279: PetscFunctionBegin;
3281: PetscAssertPointer(a, 8);
3283: dummy = (void *)(*a + mstart);
3284: PetscCall(PetscFree(dummy));
3285: PetscCall(VecRestoreArray(x, NULL));
3286: *a = NULL;
3287: PetscFunctionReturn(PETSC_SUCCESS);
3288: }
3290: /*@C
3291: VecRestoreArray3dWrite - Restores a vector after `VecGetArray3dWrite()` has been called.
3293: Logically Collective
3295: Input Parameters:
3296: + x - the vector
3297: . m - first dimension of three dimensional array
3298: . n - second dimension of the three dimensional array
3299: . p - third dimension of the three dimensional array
3300: . mstart - first index you will use in first coordinate direction (often 0)
3301: . nstart - first index in the second coordinate direction (often 0)
3302: . pstart - first index in the third coordinate direction (often 0)
3303: - a - location of pointer to array obtained from VecGetArray3d()
3305: Level: developer
3307: Notes:
3308: For regular PETSc vectors this routine does not involve any copies. For
3309: any special vectors that do not store local vector data in a contiguous
3310: array, this routine will copy the data back into the underlying
3311: vector data structure from the array obtained with `VecGetArray()`.
3313: This routine actually zeros out the `a` pointer.
3315: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3316: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3317: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3318: @*/
3319: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3320: {
3321: void *dummy;
3323: PetscFunctionBegin;
3325: PetscAssertPointer(a, 8);
3327: dummy = (void *)(*a + mstart);
3328: PetscCall(PetscFree(dummy));
3329: PetscCall(VecRestoreArrayWrite(x, NULL));
3330: *a = NULL;
3331: PetscFunctionReturn(PETSC_SUCCESS);
3332: }
3334: /*@C
3335: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3336: processor's portion of the vector data. You MUST call `VecRestoreArray4d()`
3337: when you no longer need access to the array.
3339: Logically Collective
3341: Input Parameters:
3342: + x - the vector
3343: . m - first dimension of four dimensional array
3344: . n - second dimension of four dimensional array
3345: . p - third dimension of four dimensional array
3346: . q - fourth dimension of four dimensional array
3347: . mstart - first index you will use in first coordinate direction (often 0)
3348: . nstart - first index in the second coordinate direction (often 0)
3349: . pstart - first index in the third coordinate direction (often 0)
3350: - qstart - first index in the fourth coordinate direction (often 0)
3352: Output Parameter:
3353: . a - location to put pointer to the array
3355: Level: developer
3357: Notes:
3358: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3359: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3360: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3361: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3363: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3365: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3366: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3367: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3368: @*/
3369: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3370: {
3371: PetscInt i, N, j, k;
3372: PetscScalar *aa, ***b, **c;
3374: PetscFunctionBegin;
3376: PetscAssertPointer(a, 10);
3378: PetscCall(VecGetLocalSize(x, &N));
3379: 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);
3380: PetscCall(VecGetArray(x, &aa));
3382: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3383: b = (PetscScalar ***)((*a) + m);
3384: c = (PetscScalar **)(b + m * n);
3385: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3386: for (i = 0; i < m; i++)
3387: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3388: for (i = 0; i < m; i++)
3389: for (j = 0; j < n; j++)
3390: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3391: *a -= mstart;
3392: PetscFunctionReturn(PETSC_SUCCESS);
3393: }
3395: /*@C
3396: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3397: processor's portion of the vector data. You MUST call `VecRestoreArray4dWrite()`
3398: when you no longer need access to the array.
3400: Logically Collective
3402: Input Parameters:
3403: + x - the vector
3404: . m - first dimension of four dimensional array
3405: . n - second dimension of four dimensional array
3406: . p - third dimension of four dimensional array
3407: . q - fourth dimension of four dimensional array
3408: . mstart - first index you will use in first coordinate direction (often 0)
3409: . nstart - first index in the second coordinate direction (often 0)
3410: . pstart - first index in the third coordinate direction (often 0)
3411: - qstart - first index in the fourth coordinate direction (often 0)
3413: Output Parameter:
3414: . a - location to put pointer to the array
3416: Level: developer
3418: Notes:
3419: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3420: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3421: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3422: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3424: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3426: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3427: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3428: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3429: @*/
3430: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3431: {
3432: PetscInt i, N, j, k;
3433: PetscScalar *aa, ***b, **c;
3435: PetscFunctionBegin;
3437: PetscAssertPointer(a, 10);
3439: PetscCall(VecGetLocalSize(x, &N));
3440: 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);
3441: PetscCall(VecGetArrayWrite(x, &aa));
3443: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3444: b = (PetscScalar ***)((*a) + m);
3445: c = (PetscScalar **)(b + m * n);
3446: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3447: for (i = 0; i < m; i++)
3448: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3449: for (i = 0; i < m; i++)
3450: for (j = 0; j < n; j++)
3451: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3452: *a -= mstart;
3453: PetscFunctionReturn(PETSC_SUCCESS);
3454: }
3456: /*@C
3457: VecRestoreArray4d - Restores a vector after `VecGetArray4d()` has been called.
3459: Logically Collective
3461: Input Parameters:
3462: + x - the vector
3463: . m - first dimension of four dimensional array
3464: . n - second dimension of the four dimensional array
3465: . p - third dimension of the four dimensional array
3466: . q - fourth dimension of the four dimensional array
3467: . mstart - first index you will use in first coordinate direction (often 0)
3468: . nstart - first index in the second coordinate direction (often 0)
3469: . pstart - first index in the third coordinate direction (often 0)
3470: . qstart - first index in the fourth coordinate direction (often 0)
3471: - a - location of pointer to array obtained from VecGetArray4d()
3473: Level: developer
3475: Notes:
3476: For regular PETSc vectors this routine does not involve any copies. For
3477: any special vectors that do not store local vector data in a contiguous
3478: array, this routine will copy the data back into the underlying
3479: vector data structure from the array obtained with `VecGetArray()`.
3481: This routine actually zeros out the `a` pointer.
3483: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3484: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3485: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`
3486: @*/
3487: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3488: {
3489: void *dummy;
3491: PetscFunctionBegin;
3493: PetscAssertPointer(a, 10);
3495: dummy = (void *)(*a + mstart);
3496: PetscCall(PetscFree(dummy));
3497: PetscCall(VecRestoreArray(x, NULL));
3498: *a = NULL;
3499: PetscFunctionReturn(PETSC_SUCCESS);
3500: }
3502: /*@C
3503: VecRestoreArray4dWrite - Restores a vector after `VecGetArray4dWrite()` has been called.
3505: Logically Collective
3507: Input Parameters:
3508: + x - the vector
3509: . m - first dimension of four dimensional array
3510: . n - second dimension of the four dimensional array
3511: . p - third dimension of the four dimensional array
3512: . q - fourth dimension of the four dimensional array
3513: . mstart - first index you will use in first coordinate direction (often 0)
3514: . nstart - first index in the second coordinate direction (often 0)
3515: . pstart - first index in the third coordinate direction (often 0)
3516: . qstart - first index in the fourth coordinate direction (often 0)
3517: - a - location of pointer to array obtained from `VecGetArray4d()`
3519: Level: developer
3521: Notes:
3522: For regular PETSc vectors this routine does not involve any copies. For
3523: any special vectors that do not store local vector data in a contiguous
3524: array, this routine will copy the data back into the underlying
3525: vector data structure from the array obtained with `VecGetArray()`.
3527: This routine actually zeros out the `a` pointer.
3529: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3530: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3531: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3532: @*/
3533: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3534: {
3535: void *dummy;
3537: PetscFunctionBegin;
3539: PetscAssertPointer(a, 10);
3541: dummy = (void *)(*a + mstart);
3542: PetscCall(PetscFree(dummy));
3543: PetscCall(VecRestoreArrayWrite(x, NULL));
3544: *a = NULL;
3545: PetscFunctionReturn(PETSC_SUCCESS);
3546: }
3548: /*@C
3549: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3550: processor's portion of the vector data. You MUST call `VecRestoreArray2dRead()`
3551: when you no longer need access to the array.
3553: Logically Collective
3555: Input Parameters:
3556: + x - the vector
3557: . m - first dimension of two dimensional array
3558: . n - second dimension of two dimensional array
3559: . mstart - first index you will use in first coordinate direction (often 0)
3560: - nstart - first index in the second coordinate direction (often 0)
3562: Output Parameter:
3563: . a - location to put pointer to the array
3565: Level: developer
3567: Notes:
3568: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
3569: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3570: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3571: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
3573: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3575: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3576: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3577: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3578: @*/
3579: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3580: {
3581: PetscInt i, N;
3582: const PetscScalar *aa;
3584: PetscFunctionBegin;
3586: PetscAssertPointer(a, 6);
3588: PetscCall(VecGetLocalSize(x, &N));
3589: 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);
3590: PetscCall(VecGetArrayRead(x, &aa));
3592: PetscCall(PetscMalloc1(m, a));
3593: for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3594: *a -= mstart;
3595: PetscFunctionReturn(PETSC_SUCCESS);
3596: }
3598: /*@C
3599: VecRestoreArray2dRead - Restores a vector after `VecGetArray2dRead()` has been called.
3601: Logically Collective
3603: Input Parameters:
3604: + x - the vector
3605: . m - first dimension of two dimensional array
3606: . n - second dimension of the two dimensional array
3607: . mstart - first index you will use in first coordinate direction (often 0)
3608: . nstart - first index in the second coordinate direction (often 0)
3609: - a - location of pointer to array obtained from VecGetArray2d()
3611: Level: developer
3613: Notes:
3614: For regular PETSc vectors this routine does not involve any copies. For
3615: any special vectors that do not store local vector data in a contiguous
3616: array, this routine will copy the data back into the underlying
3617: vector data structure from the array obtained with `VecGetArray()`.
3619: This routine actually zeros out the `a` pointer.
3621: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3622: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3623: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3624: @*/
3625: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3626: {
3627: void *dummy;
3629: PetscFunctionBegin;
3631: PetscAssertPointer(a, 6);
3633: dummy = (void *)(*a + mstart);
3634: PetscCall(PetscFree(dummy));
3635: PetscCall(VecRestoreArrayRead(x, NULL));
3636: *a = NULL;
3637: PetscFunctionReturn(PETSC_SUCCESS);
3638: }
3640: /*@C
3641: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3642: processor's portion of the vector data. You MUST call `VecRestoreArray1dRead()`
3643: when you no longer need access to the array.
3645: Logically Collective
3647: Input Parameters:
3648: + x - the vector
3649: . m - first dimension of two dimensional array
3650: - mstart - first index you will use in first coordinate direction (often 0)
3652: Output Parameter:
3653: . a - location to put pointer to the array
3655: Level: developer
3657: Notes:
3658: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3659: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3660: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3662: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3664: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3665: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3666: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3667: @*/
3668: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3669: {
3670: PetscInt N;
3672: PetscFunctionBegin;
3674: PetscAssertPointer(a, 4);
3676: PetscCall(VecGetLocalSize(x, &N));
3677: 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);
3678: PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3679: *a -= mstart;
3680: PetscFunctionReturn(PETSC_SUCCESS);
3681: }
3683: /*@C
3684: VecRestoreArray1dRead - Restores a vector after `VecGetArray1dRead()` has been called.
3686: Logically Collective
3688: Input Parameters:
3689: + x - the vector
3690: . m - first dimension of two dimensional array
3691: . mstart - first index you will use in first coordinate direction (often 0)
3692: - a - location of pointer to array obtained from `VecGetArray1dRead()`
3694: Level: developer
3696: Notes:
3697: For regular PETSc vectors this routine does not involve any copies. For
3698: any special vectors that do not store local vector data in a contiguous
3699: array, this routine will copy the data back into the underlying
3700: vector data structure from the array obtained with `VecGetArray1dRead()`.
3702: This routine actually zeros out the `a` pointer.
3704: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3705: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3706: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3707: @*/
3708: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3709: {
3710: PetscFunctionBegin;
3713: PetscCall(VecRestoreArrayRead(x, NULL));
3714: *a = NULL;
3715: PetscFunctionReturn(PETSC_SUCCESS);
3716: }
3718: /*@C
3719: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3720: processor's portion of the vector data. You MUST call `VecRestoreArray3dRead()`
3721: when you no longer need access to the array.
3723: Logically Collective
3725: Input Parameters:
3726: + x - the vector
3727: . m - first dimension of three dimensional array
3728: . n - second dimension of three dimensional array
3729: . p - third dimension of three dimensional array
3730: . mstart - first index you will use in first coordinate direction (often 0)
3731: . nstart - first index in the second coordinate direction (often 0)
3732: - pstart - first index in the third coordinate direction (often 0)
3734: Output Parameter:
3735: . a - location to put pointer to the array
3737: Level: developer
3739: Notes:
3740: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3741: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3742: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3743: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3dRead()`.
3745: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3747: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3748: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3749: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3750: @*/
3751: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3752: {
3753: PetscInt i, N, j;
3754: const PetscScalar *aa;
3755: PetscScalar **b;
3757: PetscFunctionBegin;
3759: PetscAssertPointer(a, 8);
3761: PetscCall(VecGetLocalSize(x, &N));
3762: 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);
3763: PetscCall(VecGetArrayRead(x, &aa));
3765: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3766: b = (PetscScalar **)((*a) + m);
3767: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3768: for (i = 0; i < m; i++)
3769: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset((PetscScalar *)aa, i * n * p + j * p - pstart);
3770: *a -= mstart;
3771: PetscFunctionReturn(PETSC_SUCCESS);
3772: }
3774: /*@C
3775: VecRestoreArray3dRead - Restores a vector after `VecGetArray3dRead()` has been called.
3777: Logically Collective
3779: Input Parameters:
3780: + x - the vector
3781: . m - first dimension of three dimensional array
3782: . n - second dimension of the three dimensional array
3783: . p - third dimension of the three dimensional array
3784: . mstart - first index you will use in first coordinate direction (often 0)
3785: . nstart - first index in the second coordinate direction (often 0)
3786: . pstart - first index in the third coordinate direction (often 0)
3787: - a - location of pointer to array obtained from `VecGetArray3dRead()`
3789: Level: developer
3791: Notes:
3792: For regular PETSc vectors this routine does not involve any copies. For
3793: any special vectors that do not store local vector data in a contiguous
3794: array, this routine will copy the data back into the underlying
3795: vector data structure from the array obtained with `VecGetArray()`.
3797: This routine actually zeros out the `a` pointer.
3799: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3800: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3801: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3802: @*/
3803: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3804: {
3805: void *dummy;
3807: PetscFunctionBegin;
3809: PetscAssertPointer(a, 8);
3811: dummy = (void *)(*a + mstart);
3812: PetscCall(PetscFree(dummy));
3813: PetscCall(VecRestoreArrayRead(x, NULL));
3814: *a = NULL;
3815: PetscFunctionReturn(PETSC_SUCCESS);
3816: }
3818: /*@C
3819: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3820: processor's portion of the vector data. You MUST call `VecRestoreArray4dRead()`
3821: when you no longer need access to the array.
3823: Logically Collective
3825: Input Parameters:
3826: + x - the vector
3827: . m - first dimension of four dimensional array
3828: . n - second dimension of four dimensional array
3829: . p - third dimension of four dimensional array
3830: . q - fourth dimension of four dimensional array
3831: . mstart - first index you will use in first coordinate direction (often 0)
3832: . nstart - first index in the second coordinate direction (often 0)
3833: . pstart - first index in the third coordinate direction (often 0)
3834: - qstart - first index in the fourth coordinate direction (often 0)
3836: Output Parameter:
3837: . a - location to put pointer to the array
3839: Level: beginner
3841: Notes:
3842: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3843: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3844: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3845: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3847: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3849: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3850: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3851: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3852: @*/
3853: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3854: {
3855: PetscInt i, N, j, k;
3856: const PetscScalar *aa;
3857: PetscScalar ***b, **c;
3859: PetscFunctionBegin;
3861: PetscAssertPointer(a, 10);
3863: PetscCall(VecGetLocalSize(x, &N));
3864: 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);
3865: PetscCall(VecGetArrayRead(x, &aa));
3867: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3868: b = (PetscScalar ***)((*a) + m);
3869: c = (PetscScalar **)(b + m * n);
3870: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3871: for (i = 0; i < m; i++)
3872: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3873: for (i = 0; i < m; i++)
3874: for (j = 0; j < n; j++)
3875: 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;
3876: *a -= mstart;
3877: PetscFunctionReturn(PETSC_SUCCESS);
3878: }
3880: /*@C
3881: VecRestoreArray4dRead - Restores a vector after `VecGetArray4d()` has been called.
3883: Logically Collective
3885: Input Parameters:
3886: + x - the vector
3887: . m - first dimension of four dimensional array
3888: . n - second dimension of the four dimensional array
3889: . p - third dimension of the four dimensional array
3890: . q - fourth dimension of the four dimensional array
3891: . mstart - first index you will use in first coordinate direction (often 0)
3892: . nstart - first index in the second coordinate direction (often 0)
3893: . pstart - first index in the third coordinate direction (often 0)
3894: . qstart - first index in the fourth coordinate direction (often 0)
3895: - a - location of pointer to array obtained from `VecGetArray4dRead()`
3897: Level: beginner
3899: Notes:
3900: For regular PETSc vectors this routine does not involve any copies. For
3901: any special vectors that do not store local vector data in a contiguous
3902: array, this routine will copy the data back into the underlying
3903: vector data structure from the array obtained with `VecGetArray()`.
3905: This routine actually zeros out the `a` pointer.
3907: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3908: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3909: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3910: @*/
3911: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3912: {
3913: void *dummy;
3915: PetscFunctionBegin;
3917: PetscAssertPointer(a, 10);
3919: dummy = (void *)(*a + mstart);
3920: PetscCall(PetscFree(dummy));
3921: PetscCall(VecRestoreArrayRead(x, NULL));
3922: *a = NULL;
3923: PetscFunctionReturn(PETSC_SUCCESS);
3924: }
3926: /*@
3927: VecLockGet - Get the current lock status of a vector
3929: Logically Collective
3931: Input Parameter:
3932: . x - the vector
3934: Output Parameter:
3935: . state - greater than zero indicates the vector is locked for read; less than zero indicates the vector is
3936: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3938: Level: advanced
3940: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3941: @*/
3942: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3943: {
3944: PetscFunctionBegin;
3946: PetscAssertPointer(state, 2);
3947: *state = x->lock;
3948: PetscFunctionReturn(PETSC_SUCCESS);
3949: }
3951: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3952: {
3953: PetscFunctionBegin;
3955: PetscAssertPointer(file, 2);
3956: PetscAssertPointer(func, 3);
3957: PetscAssertPointer(line, 4);
3958: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3959: {
3960: const int index = x->lockstack.currentsize - 1;
3962: *file = index < 0 ? NULL : x->lockstack.file[index];
3963: *func = index < 0 ? NULL : x->lockstack.function[index];
3964: *line = index < 0 ? 0 : x->lockstack.line[index];
3965: }
3966: #else
3967: *file = NULL;
3968: *func = NULL;
3969: *line = 0;
3970: #endif
3971: PetscFunctionReturn(PETSC_SUCCESS);
3972: }
3974: /*@
3975: VecLockReadPush - Push a read-only lock on a vector to prevent it from being written to
3977: Logically Collective
3979: Input Parameter:
3980: . x - the vector
3982: Level: intermediate
3984: Notes:
3985: If this is set then calls to `VecGetArray()` or `VecSetValues()` or any other routines that change the vectors values will generate an error.
3987: The call can be nested, i.e., called multiple times on the same vector, but each `VecLockReadPush()` has to have one matching
3988: `VecLockReadPop()`, which removes the latest read-only lock.
3990: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3991: @*/
3992: PetscErrorCode VecLockReadPush(Vec x)
3993: {
3994: PetscFunctionBegin;
3996: 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");
3997: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3998: {
3999: const char *file, *func;
4000: int index, line;
4002: if ((index = petscstack.currentsize - 2) < 0) {
4003: // vec was locked "outside" of petsc, either in user-land or main. the error message will
4004: // now show this function as the culprit, but it will include the stacktrace
4005: file = "unknown user-file";
4006: func = "unknown_user_function";
4007: line = 0;
4008: } else {
4009: file = petscstack.file[index];
4010: func = petscstack.function[index];
4011: line = petscstack.line[index];
4012: }
4013: PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
4014: }
4015: #endif
4016: PetscFunctionReturn(PETSC_SUCCESS);
4017: }
4019: /*@
4020: VecLockReadPop - Pop a read-only lock from a vector
4022: Logically Collective
4024: Input Parameter:
4025: . x - the vector
4027: Level: intermediate
4029: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
4030: @*/
4031: PetscErrorCode VecLockReadPop(Vec x)
4032: {
4033: PetscFunctionBegin;
4035: PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
4036: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
4037: {
4038: const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];
4040: PetscStackPop_Private(x->lockstack, previous);
4041: }
4042: #endif
4043: PetscFunctionReturn(PETSC_SUCCESS);
4044: }
4046: /*@
4047: VecLockWriteSet - Lock or unlock a vector for exclusive read/write access
4049: Logically Collective
4051: Input Parameters:
4052: + x - the vector
4053: - flg - `PETSC_TRUE` to lock the vector for exclusive read/write access; `PETSC_FALSE` to unlock it.
4055: Level: intermediate
4057: Notes:
4058: The function is useful in split-phase computations, which usually have a begin phase and an end phase.
4059: One can call `VecLockWriteSet`(x,`PETSC_TRUE`) in the begin phase to lock a vector for exclusive
4060: access, and call `VecLockWriteSet`(x,`PETSC_FALSE`) in the end phase to unlock the vector from exclusive
4061: access. In this way, one is ensured no other operations can access the vector in between. The code may like
4063: .vb
4064: VecGetArray(x,&xdata); // begin phase
4065: VecLockWriteSet(v,PETSC_TRUE);
4067: Other operations, which can not access x anymore (they can access xdata, of course)
4069: VecRestoreArray(x,&vdata); // end phase
4070: VecLockWriteSet(v,PETSC_FALSE);
4071: .ve
4073: The call can not be nested on the same vector, in other words, one can not call `VecLockWriteSet`(x,`PETSC_TRUE`)
4074: again before calling `VecLockWriteSet`(v,`PETSC_FALSE`).
4076: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
4077: @*/
4078: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
4079: {
4080: PetscFunctionBegin;
4082: if (flg) {
4083: 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");
4084: 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");
4085: x->lock = -1;
4086: } else {
4087: 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");
4088: x->lock = 0;
4089: }
4090: PetscFunctionReturn(PETSC_SUCCESS);
4091: }