Actual source code: rvector.c
1: /*
2: Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
3: These are the vector functions the user calls.
4: */
5: #include <petsc/private/vecimpl.h>
7: PetscInt VecGetSubVectorSavedStateId = -1;
9: #if PetscDefined(USE_DEBUG)
10: // this is a no-op '0' macro in optimized builds
11: PetscErrorCode VecValidValues_Internal(Vec vec, PetscInt argnum, PetscBool begin)
12: {
13: PetscFunctionBegin;
14: if (vec->petscnative || vec->ops->getarray) {
15: PetscInt n;
16: const PetscScalar *x;
17: PetscOffloadMask mask;
19: PetscCall(VecGetOffloadMask(vec, &mask));
20: if (!PetscOffloadHost(mask)) PetscFunctionReturn(PETSC_SUCCESS);
21: PetscCall(VecGetLocalSize(vec, &n));
22: PetscCall(VecGetArrayRead(vec, &x));
23: for (PetscInt i = 0; i < n; i++) {
24: if (begin) {
25: PetscCheck(!PetscIsInfOrNanScalar(x[i]), PETSC_COMM_SELF, PETSC_ERR_FP, "Vec entry at local location %" PetscInt_FMT " is not-a-number or infinite at beginning of function: Parameter number %" PetscInt_FMT, i, argnum);
26: } else {
27: PetscCheck(!PetscIsInfOrNanScalar(x[i]), PETSC_COMM_SELF, PETSC_ERR_FP, "Vec entry at local location %" PetscInt_FMT " is not-a-number or infinite at end of function: Parameter number %" PetscInt_FMT, i, argnum);
28: }
29: }
30: PetscCall(VecRestoreArrayRead(vec, &x));
31: }
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
34: #endif
36: /*@
37: VecMaxPointwiseDivide - Computes the maximum of the componentwise division `max = max_i abs(x[i]/y[i])`.
39: Logically Collective
41: Input Parameters:
42: + x - the numerators
43: - y - the denominators
45: Output Parameter:
46: . max - the result
48: Level: advanced
50: Notes:
51: `x` and `y` may be the same vector
53: if a particular `y[i]` is zero, it is treated as 1 in the above formula
55: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`
56: @*/
57: PetscErrorCode VecMaxPointwiseDivide(Vec x, Vec y, PetscReal *max)
58: {
59: PetscFunctionBegin;
62: PetscAssertPointer(max, 3);
65: PetscCheckSameTypeAndComm(x, 1, y, 2);
66: VecCheckSameSize(x, 1, y, 2);
67: VecCheckAssembled(x);
68: VecCheckAssembled(y);
69: PetscCall(VecLockReadPush(x));
70: PetscCall(VecLockReadPush(y));
71: PetscUseTypeMethod(x, maxpointwisedivide, y, max);
72: PetscCall(VecLockReadPop(x));
73: PetscCall(VecLockReadPop(y));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: /*@
78: VecDot - Computes the vector dot product.
80: Collective
82: Input Parameters:
83: + x - first vector
84: - y - second vector
86: Output Parameter:
87: . val - the dot product
89: Level: intermediate
91: Notes for Users of Complex Numbers:
92: For complex vectors, `VecDot()` computes
93: $ val = (x,y) = y^H x,
94: where y^H denotes the conjugate transpose of y. Note that this corresponds to the usual "mathematicians" complex
95: inner product where the SECOND argument gets the complex conjugate. Since the `BLASdot()` complex conjugates the first
96: first argument we call the `BLASdot()` with the arguments reversed.
98: Use `VecTDot()` for the indefinite form
99: $ val = (x,y) = y^T x,
100: where y^T denotes the transpose of y.
102: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDotRealPart()`
103: @*/
104: PetscErrorCode VecDot(Vec x, Vec y, PetscScalar *val)
105: {
106: PetscFunctionBegin;
109: PetscAssertPointer(val, 3);
112: PetscCheckSameTypeAndComm(x, 1, y, 2);
113: VecCheckSameSize(x, 1, y, 2);
114: VecCheckAssembled(x);
115: VecCheckAssembled(y);
117: PetscCall(VecLockReadPush(x));
118: PetscCall(VecLockReadPush(y));
119: PetscCall(PetscLogEventBegin(VEC_Dot, x, y, 0, 0));
120: PetscUseTypeMethod(x, dot, y, val);
121: PetscCall(PetscLogEventEnd(VEC_Dot, x, y, 0, 0));
122: PetscCall(VecLockReadPop(x));
123: PetscCall(VecLockReadPop(y));
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: /*@
128: VecDotRealPart - Computes the real part of the vector dot product.
130: Collective
132: Input Parameters:
133: + x - first vector
134: - y - second vector
136: Output Parameter:
137: . val - the real part of the dot product;
139: Level: intermediate
141: Notes for Users of Complex Numbers:
142: See `VecDot()` for more details on the definition of the dot product for complex numbers
144: For real numbers this returns the same value as `VecDot()`
146: For complex numbers in C^n (that is a vector of n components with a complex number for each component) this is equal to the usual real dot product on the
147: the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)
149: Developer Notes:
150: This is not currently optimized to compute only the real part of the dot product.
152: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDot()`, `VecDotNorm2()`
153: @*/
154: PetscErrorCode VecDotRealPart(Vec x, Vec y, PetscReal *val)
155: {
156: PetscScalar fdot;
158: PetscFunctionBegin;
159: PetscCall(VecDot(x, y, &fdot));
160: *val = PetscRealPart(fdot);
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: /*@
165: VecNorm - Computes the vector norm.
167: Collective
169: Input Parameters:
170: + x - the vector
171: - type - the type of the norm requested
173: Output Parameter:
174: . val - the norm
176: Level: intermediate
178: Notes:
179: See `NormType` for descriptions of each norm.
181: For complex numbers `NORM_1` will return the traditional 1 norm of the 2 norm of the complex
182: numbers; that is the 1 norm of the absolute values of the complex entries. In PETSc 3.6 and
183: earlier releases it returned the 1 norm of the 1 norm of the complex entries (what is
184: returned by the BLAS routine `asum()`). Both are valid norms but most people expect the former.
186: This routine stashes the computed norm value, repeated calls before the vector entries are
187: changed are then rapid since the precomputed value is immediately available. Certain vector
188: operations such as `VecSet()` store the norms so the value is immediately available and does
189: not need to be explicitly computed. `VecScale()` updates any stashed norm values, thus calls
190: after `VecScale()` do not need to explicitly recompute the norm.
192: .seealso: [](ch_vectors), `Vec`, `NormType`, `VecDot()`, `VecTDot()`, `VecDotBegin()`, `VecDotEnd()`, `VecNormAvailable()`,
193: `VecNormBegin()`, `VecNormEnd()`, `NormType()`
194: @*/
195: PetscErrorCode VecNorm(Vec x, NormType type, PetscReal *val)
196: {
197: PetscBool flg = PETSC_TRUE;
199: PetscFunctionBegin;
202: VecCheckAssembled(x);
204: PetscAssertPointer(val, 3);
206: PetscCall(VecNormAvailable(x, type, &flg, val));
207: // check that all MPI processes call this routine together and have same availability
208: if (PetscDefined(USE_DEBUG)) {
209: PetscMPIInt b0 = (PetscMPIInt)flg, b1[2], b2[2];
210: b1[0] = -b0;
211: b1[1] = b0;
212: PetscCall(MPIU_Allreduce(b1, b2, 2, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
213: PetscCheck(-b2[0] == b2[1], PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_WRONGSTATE, "Some MPI processes have cached %s norm, others do not. This may happen when some MPI processes call VecGetArray() and some others do not.", NormTypes[type]);
214: if (flg) {
215: PetscReal b1[2], b2[2];
216: b1[0] = -(*val);
217: b1[1] = *val;
218: PetscCall(MPIU_Allreduce(b1, b2, 2, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)x)));
219: PetscCheck(-b2[0] == b2[1], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Difference in cached %s norms: local %g", NormTypes[type], (double)*val);
220: }
221: }
222: if (flg) PetscFunctionReturn(PETSC_SUCCESS);
224: PetscCall(VecLockReadPush(x));
225: PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
226: PetscUseTypeMethod(x, norm, type, val);
227: PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
228: PetscCall(VecLockReadPop(x));
230: if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: /*@
235: VecNormAvailable - Returns the vector norm if it is already known. That is, it has been previously computed and cached in the vector
237: Not Collective
239: Input Parameters:
240: + x - the vector
241: - type - one of `NORM_1` (sum_i |x[i]|), `NORM_2` sqrt(sum_i (x[i])^2), `NORM_INFINITY` max_i |x[i]|. Also available
242: `NORM_1_AND_2`, which computes both norms and stores them
243: in a two element array.
245: Output Parameters:
246: + available - `PETSC_TRUE` if the val returned is valid
247: - val - the norm
249: Level: intermediate
251: Developer Notes:
252: `PETSC_HAVE_SLOW_BLAS_NORM2` will cause a C (loop unrolled) version of the norm to be used, rather
253: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
254: (as opposed to vendor provided) because the FORTRAN BLAS `NRM2()` routine is very slow.
256: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
257: `VecNormBegin()`, `VecNormEnd()`
258: @*/
259: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
260: {
261: PetscFunctionBegin;
264: PetscAssertPointer(available, 3);
265: PetscAssertPointer(val, 4);
267: if (type == NORM_1_AND_2) {
268: *available = PETSC_FALSE;
269: } else {
270: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
271: }
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: /*@
276: VecNormalize - Normalizes a vector by its 2-norm.
278: Collective
280: Input Parameter:
281: . x - the vector
283: Output Parameter:
284: . val - the vector norm before normalization. May be `NULL` if the value is not needed.
286: Level: intermediate
288: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
289: @*/
290: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
291: {
292: PetscReal norm;
294: PetscFunctionBegin;
297: PetscCall(VecSetErrorIfLocked(x, 1));
298: if (val) PetscAssertPointer(val, 2);
299: PetscCall(PetscLogEventBegin(VEC_Normalize, x, 0, 0, 0));
300: PetscCall(VecNorm(x, NORM_2, &norm));
301: if (norm == 0.0) PetscCall(PetscInfo(x, "Vector of zero norm can not be normalized; Returning only the zero norm\n"));
302: else if (PetscIsInfOrNanReal(norm)) PetscCall(PetscInfo(x, "Vector with Inf or Nan norm can not be normalized; Returning only the norm\n"));
303: else {
304: PetscScalar s = 1.0 / norm;
305: PetscCall(VecScale(x, s));
306: }
307: PetscCall(PetscLogEventEnd(VEC_Normalize, x, 0, 0, 0));
308: if (val) *val = norm;
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: /*@C
313: VecMax - Determines the vector component with maximum real part and its location.
315: Collective
317: Input Parameter:
318: . x - the vector
320: Output Parameters:
321: + p - the index of `val` (pass `NULL` if you don't want this) in the vector
322: - val - the maximum component
324: Level: intermediate
326: Notes:
327: Returns the value `PETSC_MIN_REAL` and negative `p` if the vector is of length 0.
329: Returns the smallest index with the maximum value
331: Developer Note:
332: The Nag Fortran compiler does not like the symbol name VecMax
334: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
335: @*/
336: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
337: {
338: PetscFunctionBegin;
341: VecCheckAssembled(x);
342: if (p) PetscAssertPointer(p, 2);
343: PetscAssertPointer(val, 3);
344: PetscCall(VecLockReadPush(x));
345: PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
346: PetscUseTypeMethod(x, max, p, val);
347: PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
348: PetscCall(VecLockReadPop(x));
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: /*@C
353: VecMin - Determines the vector component with minimum real part and its location.
355: Collective
357: Input Parameter:
358: . x - the vector
360: Output Parameters:
361: + p - the index of `val` (pass `NULL` if you don't want this location) in the vector
362: - val - the minimum component
364: Level: intermediate
366: Notes:
367: Returns the value `PETSC_MAX_REAL` and negative `p` if the vector is of length 0.
369: This returns the smallest index with the minimum value
371: Developer Note:
372: The Nag Fortran compiler does not like the symbol name VecMin
374: .seealso: [](ch_vectors), `Vec`, `VecMax()`
375: @*/
376: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
377: {
378: PetscFunctionBegin;
381: VecCheckAssembled(x);
382: if (p) PetscAssertPointer(p, 2);
383: PetscAssertPointer(val, 3);
384: PetscCall(VecLockReadPush(x));
385: PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
386: PetscUseTypeMethod(x, min, p, val);
387: PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
388: PetscCall(VecLockReadPop(x));
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: /*@
393: VecTDot - Computes an indefinite vector dot product. That is, this
394: routine does NOT use the complex conjugate.
396: Collective
398: Input Parameters:
399: + x - first vector
400: - y - second vector
402: Output Parameter:
403: . val - the dot product
405: Level: intermediate
407: Notes for Users of Complex Numbers:
408: For complex vectors, `VecTDot()` computes the indefinite form
409: $ val = (x,y) = y^T x,
410: where y^T denotes the transpose of y.
412: Use `VecDot()` for the inner product
413: $ val = (x,y) = y^H x,
414: where y^H denotes the conjugate transpose of y.
416: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
417: @*/
418: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
419: {
420: PetscFunctionBegin;
423: PetscAssertPointer(val, 3);
426: PetscCheckSameTypeAndComm(x, 1, y, 2);
427: VecCheckSameSize(x, 1, y, 2);
428: VecCheckAssembled(x);
429: VecCheckAssembled(y);
431: PetscCall(VecLockReadPush(x));
432: PetscCall(VecLockReadPush(y));
433: PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
434: PetscUseTypeMethod(x, tdot, y, val);
435: PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
436: PetscCall(VecLockReadPop(x));
437: PetscCall(VecLockReadPop(y));
438: PetscFunctionReturn(PETSC_SUCCESS);
439: }
441: PetscErrorCode VecScaleAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
442: {
443: PetscReal norms[4];
444: PetscBool flgs[4];
445: PetscScalar one = 1.0;
447: PetscFunctionBegin;
450: VecCheckAssembled(x);
451: PetscCall(VecSetErrorIfLocked(x, 1));
453: if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);
455: /* get current stashed norms */
456: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
458: PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
459: VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
460: PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));
462: PetscCall(PetscObjectStateIncrease((PetscObject)x));
463: /* put the scaled stashed norms back into the Vec */
464: for (PetscInt i = 0; i < 4; i++) {
465: PetscReal ar = PetscAbsScalar(alpha);
466: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
467: }
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
471: /*@
472: VecScale - Scales a vector.
474: Logically Collective
476: Input Parameters:
477: + x - the vector
478: - alpha - the scalar
480: Level: intermediate
482: Note:
483: For a vector with n components, `VecScale()` computes x[i] = alpha * x[i], for i=1,...,n.
485: .seealso: [](ch_vectors), `Vec`, `VecSet()`
486: @*/
487: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
488: {
489: PetscFunctionBegin;
490: PetscCall(VecScaleAsync_Private(x, alpha, NULL));
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
495: {
496: PetscFunctionBegin;
499: VecCheckAssembled(x);
501: PetscCall(VecSetErrorIfLocked(x, 1));
503: if (alpha == 0) {
504: PetscReal norm;
505: PetscBool set;
507: PetscCall(VecNormAvailable(x, NORM_2, &set, &norm));
508: if (set == PETSC_TRUE && norm == 0) PetscFunctionReturn(PETSC_SUCCESS);
509: }
510: PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
511: VecMethodDispatch(x, dctx, VecAsyncFnName(Set), set, (Vec, PetscScalar, PetscDeviceContext), alpha);
512: PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
513: PetscCall(PetscObjectStateIncrease((PetscObject)x));
515: /* norms can be simply set (if |alpha|*N not too large) */
516: {
517: PetscReal val = PetscAbsScalar(alpha);
518: const PetscInt N = x->map->N;
520: if (N == 0) {
521: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0));
522: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
523: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
524: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
525: } else if (val > PETSC_MAX_REAL / N) {
526: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
527: } else {
528: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
529: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
530: val *= PetscSqrtReal((PetscReal)N);
531: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
532: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
533: }
534: }
535: PetscFunctionReturn(PETSC_SUCCESS);
536: }
538: /*@
539: VecSet - Sets all components of a vector to a single scalar value.
541: Logically Collective
543: Input Parameters:
544: + x - the vector
545: - alpha - the scalar
547: Level: beginner
549: Notes:
550: For a vector of dimension n, `VecSet()` sets x[i] = alpha, for i=1,...,n,
551: so that all vector entries then equal the identical
552: scalar value, `alpha`. Use the more general routine
553: `VecSetValues()` to set different vector entries.
555: You CANNOT call this after you have called `VecSetValues()` but before you call
556: `VecAssemblyBegin()`
558: If `alpha` is zero and the norm of the vector is known to be zero then this skips the unneeded zeroing process
560: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
561: @*/
562: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
563: {
564: PetscFunctionBegin;
565: PetscCall(VecSetAsync_Private(x, alpha, NULL));
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: PetscErrorCode VecAXPYAsync_Private(Vec y, PetscScalar alpha, Vec x, PetscDeviceContext dctx)
570: {
571: PetscFunctionBegin;
576: PetscCheckSameTypeAndComm(x, 3, y, 1);
577: VecCheckSameSize(x, 3, y, 1);
578: VecCheckAssembled(x);
579: VecCheckAssembled(y);
581: if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
582: PetscCall(VecSetErrorIfLocked(y, 1));
583: if (x == y) {
584: PetscCall(VecScale(y, alpha + 1.0));
585: PetscFunctionReturn(PETSC_SUCCESS);
586: }
587: PetscCall(VecLockReadPush(x));
588: PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
589: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPY), axpy, (Vec, PetscScalar, Vec, PetscDeviceContext), alpha, x);
590: PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
591: PetscCall(VecLockReadPop(x));
592: PetscCall(PetscObjectStateIncrease((PetscObject)y));
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
595: /*@
596: VecAXPY - Computes `y = alpha x + y`.
598: Logically Collective
600: Input Parameters:
601: + alpha - the scalar
602: . x - vector scale by `alpha`
603: - y - vector accumulated into
605: Output Parameter:
606: . y - output vector
608: Level: intermediate
610: Notes:
611: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
612: .vb
613: VecAXPY(y,alpha,x) y = alpha x + y
614: VecAYPX(y,beta,x) y = x + beta y
615: VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
616: VecWAXPY(w,alpha,x,y) w = alpha x + y
617: VecAXPBYPCZ(z,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
618: VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
619: .ve
621: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
622: @*/
623: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
624: {
625: PetscFunctionBegin;
626: PetscCall(VecAXPYAsync_Private(y, alpha, x, NULL));
627: PetscFunctionReturn(PETSC_SUCCESS);
628: }
630: PetscErrorCode VecAYPXAsync_Private(Vec y, PetscScalar beta, Vec x, PetscDeviceContext dctx)
631: {
632: PetscFunctionBegin;
637: PetscCheckSameTypeAndComm(x, 3, y, 1);
638: VecCheckSameSize(x, 1, y, 3);
639: VecCheckAssembled(x);
640: VecCheckAssembled(y);
642: PetscCall(VecSetErrorIfLocked(y, 1));
643: if (x == y) {
644: PetscCall(VecScale(y, beta + 1.0));
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
647: PetscCall(VecLockReadPush(x));
648: if (beta == (PetscScalar)0.0) {
649: PetscCall(VecCopy(x, y));
650: } else {
651: PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
652: VecMethodDispatch(y, dctx, VecAsyncFnName(AYPX), aypx, (Vec, PetscScalar, Vec, PetscDeviceContext), beta, x);
653: PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
654: PetscCall(PetscObjectStateIncrease((PetscObject)y));
655: }
656: PetscCall(VecLockReadPop(x));
657: PetscFunctionReturn(PETSC_SUCCESS);
658: }
660: /*@
661: VecAYPX - Computes `y = x + beta y`.
663: Logically Collective
665: Input Parameters:
666: + beta - the scalar
667: . x - the unscaled vector
668: - y - the vector to be scaled
670: Output Parameter:
671: . y - output vector
673: Level: intermediate
675: Developer Notes:
676: The implementation is optimized for `beta` of -1.0, 0.0, and 1.0
678: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
679: @*/
680: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
681: {
682: PetscFunctionBegin;
683: PetscCall(VecAYPXAsync_Private(y, beta, x, NULL));
684: PetscFunctionReturn(PETSC_SUCCESS);
685: }
687: PetscErrorCode VecAXPBYAsync_Private(Vec y, PetscScalar alpha, PetscScalar beta, Vec x, PetscDeviceContext dctx)
688: {
689: PetscFunctionBegin;
694: PetscCheckSameTypeAndComm(x, 4, y, 1);
695: VecCheckSameSize(y, 1, x, 4);
696: VecCheckAssembled(x);
697: VecCheckAssembled(y);
700: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
701: if (x == y) {
702: PetscCall(VecScale(y, alpha + beta));
703: PetscFunctionReturn(PETSC_SUCCESS);
704: }
706: PetscCall(VecSetErrorIfLocked(y, 1));
707: PetscCall(VecLockReadPush(x));
708: PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
709: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPBY), axpby, (Vec, PetscScalar, PetscScalar, Vec, PetscDeviceContext), alpha, beta, x);
710: PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
711: PetscCall(PetscObjectStateIncrease((PetscObject)y));
712: PetscCall(VecLockReadPop(x));
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }
716: /*@
717: VecAXPBY - Computes `y = alpha x + beta y`.
719: Logically Collective
721: Input Parameters:
722: + alpha - first scalar
723: . beta - second scalar
724: . x - the first scaled vector
725: - y - the second scaled vector
727: Output Parameter:
728: . y - output vector
730: Level: intermediate
732: Developer Notes:
733: The implementation is optimized for `alpha` and/or `beta` values of 0.0 and 1.0
735: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
736: @*/
737: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
738: {
739: PetscFunctionBegin;
740: PetscCall(VecAXPBYAsync_Private(y, alpha, beta, x, NULL));
741: PetscFunctionReturn(PETSC_SUCCESS);
742: }
744: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
745: {
746: PetscFunctionBegin;
753: PetscCheckSameTypeAndComm(x, 5, y, 6);
754: PetscCheckSameTypeAndComm(x, 5, z, 1);
755: VecCheckSameSize(x, 5, y, 6);
756: VecCheckSameSize(x, 5, z, 1);
757: PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
758: PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
759: VecCheckAssembled(x);
760: VecCheckAssembled(y);
761: VecCheckAssembled(z);
765: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
767: PetscCall(VecSetErrorIfLocked(z, 1));
768: PetscCall(VecLockReadPush(x));
769: PetscCall(VecLockReadPush(y));
770: PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
771: VecMethodDispatch(z, dctx, VecAsyncFnName(AXPBYPCZ), axpbypcz, (Vec, PetscScalar, PetscScalar, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, beta, gamma, x, y);
772: PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
773: PetscCall(PetscObjectStateIncrease((PetscObject)z));
774: PetscCall(VecLockReadPop(x));
775: PetscCall(VecLockReadPop(y));
776: PetscFunctionReturn(PETSC_SUCCESS);
777: }
778: /*@
779: VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`
781: Logically Collective
783: Input Parameters:
784: + alpha - first scalar
785: . beta - second scalar
786: . gamma - third scalar
787: . x - first vector
788: . y - second vector
789: - z - third vector
791: Output Parameter:
792: . z - output vector
794: Level: intermediate
796: Note:
797: `x`, `y` and `z` must be different vectors
799: Developer Notes:
800: The implementation is optimized for `alpha` of 1.0 and `gamma` of 1.0 or 0.0
802: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
803: @*/
804: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
805: {
806: PetscFunctionBegin;
807: PetscCall(VecAXPBYPCZAsync_Private(z, alpha, beta, gamma, x, y, NULL));
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
812: {
813: PetscFunctionBegin;
820: PetscCheckSameTypeAndComm(x, 3, y, 4);
821: PetscCheckSameTypeAndComm(y, 4, w, 1);
822: VecCheckSameSize(x, 3, y, 4);
823: VecCheckSameSize(x, 3, w, 1);
824: PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
825: PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
826: VecCheckAssembled(x);
827: VecCheckAssembled(y);
829: PetscCall(VecSetErrorIfLocked(w, 1));
831: PetscCall(VecLockReadPush(x));
832: PetscCall(VecLockReadPush(y));
833: if (alpha == (PetscScalar)0.0) {
834: PetscCall(VecCopyAsync_Private(y, w, dctx));
835: } else {
836: PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
837: VecMethodDispatch(w, dctx, VecAsyncFnName(WAXPY), waxpy, (Vec, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, x, y);
838: PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
839: PetscCall(PetscObjectStateIncrease((PetscObject)w));
840: }
841: PetscCall(VecLockReadPop(x));
842: PetscCall(VecLockReadPop(y));
843: PetscFunctionReturn(PETSC_SUCCESS);
844: }
846: /*@
847: VecWAXPY - Computes `w = alpha x + y`.
849: Logically Collective
851: Input Parameters:
852: + alpha - the scalar
853: . x - first vector, multiplied by `alpha`
854: - y - second vector
856: Output Parameter:
857: . w - the result
859: Level: intermediate
861: Note:
862: `w` cannot be either `x` or `y`, but `x` and `y` can be the same
864: Developer Notes:
865: The implementation is optimized for alpha of -1.0, 0.0, and 1.0
867: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
868: @*/
869: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
870: {
871: PetscFunctionBegin;
872: PetscCall(VecWAXPYAsync_Private(w, alpha, x, y, NULL));
873: PetscFunctionReturn(PETSC_SUCCESS);
874: }
876: /*@
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 `VecSetValues()` 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 (x->map->mapping) {
1086: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1087: PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1088: PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1089: if (ni > 128) PetscCall(PetscFree(lix));
1090: } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1091: } else PetscUseTypeMethod(x, setvalueslocal, ni, ix, y, iora);
1092: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1093: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1094: PetscFunctionReturn(PETSC_SUCCESS);
1095: }
1097: /*@
1098: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1099: using a local ordering of the nodes.
1101: Not Collective
1103: Input Parameters:
1104: + x - vector to insert in
1105: . ni - number of blocks to add
1106: . ix - indices where to add in block count, not element count
1107: . y - array of values
1108: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1110: Level: intermediate
1112: Notes:
1113: `VecSetValuesBlockedLocal()` sets x[bs*ix[i]+j] = y[bs*i+j],
1114: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with `VecSetBlockSize()`.
1116: Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1117: options cannot be mixed without intervening calls to the assembly
1118: routines.
1120: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1121: MUST be called after all calls to `VecSetValuesBlockedLocal()` have been completed.
1123: `VecSetValuesBlockedLocal()` uses 0-based indices in Fortran as well as in C.
1125: Fortran Note:
1126: If any of `ix` and `y` are scalars pass them using, for example,
1127: .vb
1128: VecSetValuesBlockedLocal(mat, one, [ix], [y], INSERT_VALUES)
1129: .ve
1131: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1132: `VecSetLocalToGlobalMapping()`
1133: @*/
1134: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1135: {
1136: PetscInt lixp[128], *lix = lixp;
1138: PetscFunctionBeginHot;
1140: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1141: PetscAssertPointer(ix, 3);
1142: PetscAssertPointer(y, 4);
1144: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1145: if (x->map->mapping) {
1146: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1147: PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1148: PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1149: if (ni > 128) PetscCall(PetscFree(lix));
1150: } else {
1151: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1152: }
1153: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1154: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1155: PetscFunctionReturn(PETSC_SUCCESS);
1156: }
1158: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1159: {
1160: PetscFunctionBegin;
1163: VecCheckAssembled(x);
1165: if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1166: PetscAssertPointer(y, 3);
1167: for (PetscInt i = 0; i < nv; ++i) {
1170: PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1171: VecCheckSameSize(x, 1, y[i], 3);
1172: VecCheckAssembled(y[i]);
1173: PetscCall(VecLockReadPush(y[i]));
1174: }
1175: PetscAssertPointer(result, 4);
1178: PetscCall(VecLockReadPush(x));
1179: PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1180: PetscCall((*mxdot)(x, nv, y, result));
1181: PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1182: PetscCall(VecLockReadPop(x));
1183: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1184: PetscFunctionReturn(PETSC_SUCCESS);
1185: }
1187: /*@
1188: VecMTDot - Computes indefinite vector multiple dot products.
1189: That is, it does NOT use the complex conjugate.
1191: Collective
1193: Input Parameters:
1194: + x - one vector
1195: . nv - number of vectors
1196: - y - array of vectors. Note that vectors are pointers
1198: Output Parameter:
1199: . val - array of the dot products
1201: Level: intermediate
1203: Notes for Users of Complex Numbers:
1204: For complex vectors, `VecMTDot()` computes the indefinite form
1205: $ val = (x,y) = y^T x,
1206: where y^T denotes the transpose of y.
1208: Use `VecMDot()` for the inner product
1209: $ val = (x,y) = y^H x,
1210: where y^H denotes the conjugate transpose of y.
1212: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1213: @*/
1214: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1215: {
1216: PetscFunctionBegin;
1218: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1219: PetscFunctionReturn(PETSC_SUCCESS);
1220: }
1222: /*@
1223: VecMDot - Computes multiple vector dot products.
1225: Collective
1227: Input Parameters:
1228: + x - one vector
1229: . nv - number of vectors
1230: - y - array of vectors.
1232: Output Parameter:
1233: . val - array of the dot products (does not allocate the array)
1235: Level: intermediate
1237: Notes for Users of Complex Numbers:
1238: For complex vectors, `VecMDot()` computes
1239: $ val = (x,y) = y^H x,
1240: where y^H denotes the conjugate transpose of y.
1242: Use `VecMTDot()` for the indefinite form
1243: $ val = (x,y) = y^T x,
1244: where y^T denotes the transpose of y.
1246: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`
1247: @*/
1248: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1249: {
1250: PetscFunctionBegin;
1252: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1253: PetscFunctionReturn(PETSC_SUCCESS);
1254: }
1256: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
1257: {
1258: PetscFunctionBegin;
1260: VecCheckAssembled(y);
1262: PetscCall(VecSetErrorIfLocked(y, 1));
1263: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1264: if (nv) {
1265: PetscInt zeros = 0;
1267: PetscAssertPointer(alpha, 3);
1268: PetscAssertPointer(x, 4);
1269: for (PetscInt i = 0; i < nv; ++i) {
1273: PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1274: VecCheckSameSize(y, 1, x[i], 4);
1275: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1276: VecCheckAssembled(x[i]);
1277: PetscCall(VecLockReadPush(x[i]));
1278: zeros += alpha[i] == (PetscScalar)0.0;
1279: }
1281: if (zeros < nv) {
1282: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1283: VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1284: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1285: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1286: }
1288: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1289: }
1290: PetscFunctionReturn(PETSC_SUCCESS);
1291: }
1293: /*@
1294: VecMAXPY - Computes `y = y + sum alpha[i] x[i]`
1296: Logically Collective
1298: Input Parameters:
1299: + nv - number of scalars and x-vectors
1300: . alpha - array of scalars
1301: . y - one vector
1302: - x - array of vectors
1304: Level: intermediate
1306: Note:
1307: `y` cannot be any of the `x` vectors
1309: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`,`VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1310: @*/
1311: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1312: {
1313: PetscFunctionBegin;
1314: PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1315: PetscFunctionReturn(PETSC_SUCCESS);
1316: }
1318: /*@
1319: VecMAXPBY - Computes `y = beta y + sum alpha[i] x[i]`
1321: Logically Collective
1323: Input Parameters:
1324: + nv - number of scalars and x-vectors
1325: . alpha - array of scalars
1326: . beta - scalar
1327: . y - one vector
1328: - x - array of vectors
1330: Level: intermediate
1332: Note:
1333: `y` cannot be any of the `x` vectors.
1335: Developer Notes:
1336: This is a convenience routine, but implementations might be able to optimize it, for example, when `beta` is zero.
1338: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1339: @*/
1340: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
1341: {
1342: PetscFunctionBegin;
1344: VecCheckAssembled(y);
1346: PetscCall(VecSetErrorIfLocked(y, 1));
1347: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1350: if (y->ops->maxpby) {
1351: PetscInt zeros = 0;
1353: if (nv) {
1354: PetscAssertPointer(alpha, 3);
1355: PetscAssertPointer(x, 5);
1356: }
1358: for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
1362: PetscCheckSameTypeAndComm(y, 1, x[i], 5);
1363: VecCheckSameSize(y, 1, x[i], 5);
1364: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1365: VecCheckAssembled(x[i]);
1366: PetscCall(VecLockReadPush(x[i]));
1367: zeros += alpha[i] == (PetscScalar)0.0;
1368: }
1370: if (zeros < nv) { // has nonzero alpha
1371: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1372: PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1373: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1374: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1375: } else {
1376: PetscCall(VecScale(y, beta));
1377: }
1379: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1380: } else { // no maxpby
1381: if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1382: else PetscCall(VecScale(y, beta));
1383: PetscCall(VecMAXPY(y, nv, alpha, x));
1384: }
1385: PetscFunctionReturn(PETSC_SUCCESS);
1386: }
1388: /*@
1389: VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
1390: in the order they appear in the array. The concatenated vector resides on the same
1391: communicator and is the same type as the source vectors.
1393: Collective
1395: Input Parameters:
1396: + nx - number of vectors to be concatenated
1397: - X - array containing the vectors to be concatenated in the order of concatenation
1399: Output Parameters:
1400: + Y - concatenated vector
1401: - x_is - array of index sets corresponding to the concatenated components of `Y` (pass `NULL` if not needed)
1403: Level: advanced
1405: Notes:
1406: Concatenation is similar to the functionality of a `VECNEST` object; they both represent combination of
1407: different vector spaces. However, concatenated vectors do not store any information about their
1408: sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1409: manipulation of data in the concatenated vector that corresponds to the original components at creation.
1411: This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1412: has to operate on combined vector spaces and cannot utilize `VECNEST` objects due to incompatibility with
1413: bound projections.
1415: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1416: @*/
1417: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1418: {
1419: MPI_Comm comm;
1420: VecType vec_type;
1421: Vec Ytmp, Xtmp;
1422: IS *is_tmp;
1423: PetscInt i, shift = 0, Xnl, Xng, Xbegin;
1425: PetscFunctionBegin;
1429: PetscAssertPointer(Y, 3);
1431: if ((*X)->ops->concatenate) {
1432: /* use the dedicated concatenation function if available */
1433: PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1434: } else {
1435: /* loop over vectors and start creating IS */
1436: comm = PetscObjectComm((PetscObject)*X);
1437: PetscCall(VecGetType(*X, &vec_type));
1438: PetscCall(PetscMalloc1(nx, &is_tmp));
1439: for (i = 0; i < nx; i++) {
1440: PetscCall(VecGetSize(X[i], &Xng));
1441: PetscCall(VecGetLocalSize(X[i], &Xnl));
1442: PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1443: PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1444: shift += Xng;
1445: }
1446: /* create the concatenated vector */
1447: PetscCall(VecCreate(comm, &Ytmp));
1448: PetscCall(VecSetType(Ytmp, vec_type));
1449: PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1450: PetscCall(VecSetUp(Ytmp));
1451: /* copy data from X array to Y and return */
1452: for (i = 0; i < nx; i++) {
1453: PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1454: PetscCall(VecCopy(X[i], Xtmp));
1455: PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1456: }
1457: *Y = Ytmp;
1458: if (x_is) {
1459: *x_is = is_tmp;
1460: } else {
1461: for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1462: PetscCall(PetscFree(is_tmp));
1463: }
1464: }
1465: PetscFunctionReturn(PETSC_SUCCESS);
1466: }
1468: /* A helper function for VecGetSubVector to check if we can implement it with no-copy (i.e. the subvector shares
1469: memory with the original vector), and the block size of the subvector.
1471: Input Parameters:
1472: + X - the original vector
1473: - is - the index set of the subvector
1475: Output Parameters:
1476: + contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
1477: . start - start of contiguous block, as an offset from the start of the ownership range of the original vector
1478: - blocksize - the block size of the subvector
1480: */
1481: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1482: {
1483: PetscInt gstart, gend, lstart;
1484: PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1485: PetscInt n, N, ibs, vbs, bs = -1;
1487: PetscFunctionBegin;
1488: PetscCall(ISGetLocalSize(is, &n));
1489: PetscCall(ISGetSize(is, &N));
1490: PetscCall(ISGetBlockSize(is, &ibs));
1491: PetscCall(VecGetBlockSize(X, &vbs));
1492: PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1493: PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1494: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1495: if (ibs > 1) {
1496: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1497: bs = ibs;
1498: } else {
1499: if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1500: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1501: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1502: }
1504: *contig = red[0];
1505: *start = lstart;
1506: *blocksize = bs;
1507: PetscFunctionReturn(PETSC_SUCCESS);
1508: }
1510: /* A helper function for VecGetSubVector, to be used when we have to build a standalone subvector through VecScatter
1512: Input Parameters:
1513: + X - the original vector
1514: . is - the index set of the subvector
1515: - bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()
1517: Output Parameter:
1518: . Z - the subvector, which will compose the VecScatter context on output
1519: */
1520: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1521: {
1522: PetscInt n, N;
1523: VecScatter vscat;
1524: Vec Y;
1526: PetscFunctionBegin;
1527: PetscCall(ISGetLocalSize(is, &n));
1528: PetscCall(ISGetSize(is, &N));
1529: PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1530: PetscCall(VecSetSizes(Y, n, N));
1531: PetscCall(VecSetBlockSize(Y, bs));
1532: PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1533: PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1534: PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1535: PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1536: PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1537: PetscCall(VecScatterDestroy(&vscat));
1538: *Z = Y;
1539: PetscFunctionReturn(PETSC_SUCCESS);
1540: }
1542: /*@
1543: VecGetSubVector - Gets a vector representing part of another vector
1545: Collective
1547: Input Parameters:
1548: + X - vector from which to extract a subvector
1549: - is - index set representing portion of `X` to extract
1551: Output Parameter:
1552: . Y - subvector corresponding to `is`
1554: Level: advanced
1556: Notes:
1557: The subvector `Y` should be returned with `VecRestoreSubVector()`.
1558: `X` and `is` must be defined on the same communicator
1560: Changes to the subvector will be reflected in the `X` vector on the call to `VecRestoreSubVector()`.
1562: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1563: modifying the subvector. Other non-overlapping subvectors can still be obtained from `X` using this function.
1565: 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`.
1567: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1568: @*/
1569: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1570: {
1571: Vec Z;
1573: PetscFunctionBegin;
1576: PetscCheckSameComm(X, 1, is, 2);
1577: PetscAssertPointer(Y, 3);
1578: if (X->ops->getsubvector) {
1579: PetscUseTypeMethod(X, getsubvector, is, &Z);
1580: } else { /* Default implementation currently does no caching */
1581: PetscBool contig;
1582: PetscInt n, N, start, bs;
1584: PetscCall(ISGetLocalSize(is, &n));
1585: PetscCall(ISGetSize(is, &N));
1586: PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1587: if (contig) { /* We can do a no-copy implementation */
1588: const PetscScalar *x;
1589: PetscInt state = 0;
1590: PetscBool isstd, iscuda, iship;
1592: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1593: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1594: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1595: if (iscuda) {
1596: #if defined(PETSC_HAVE_CUDA)
1597: const PetscScalar *x_d;
1598: PetscMPIInt size;
1599: PetscOffloadMask flg;
1601: PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1602: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1603: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1604: if (x) x += start;
1605: if (x_d) x_d += start;
1606: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1607: if (size == 1) {
1608: PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1609: } else {
1610: PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1611: }
1612: Z->offloadmask = flg;
1613: #endif
1614: } else if (iship) {
1615: #if defined(PETSC_HAVE_HIP)
1616: const PetscScalar *x_d;
1617: PetscMPIInt size;
1618: PetscOffloadMask flg;
1620: PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1621: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1622: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1623: if (x) x += start;
1624: if (x_d) x_d += start;
1625: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1626: if (size == 1) {
1627: PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1628: } else {
1629: PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1630: }
1631: Z->offloadmask = flg;
1632: #endif
1633: } else if (isstd) {
1634: PetscMPIInt size;
1636: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1637: PetscCall(VecGetArrayRead(X, &x));
1638: if (x) x += start;
1639: if (size == 1) {
1640: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1641: } else {
1642: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1643: }
1644: PetscCall(VecRestoreArrayRead(X, &x));
1645: } else { /* default implementation: use place array */
1646: PetscCall(VecGetArrayRead(X, &x));
1647: PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1648: PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1649: PetscCall(VecSetSizes(Z, n, N));
1650: PetscCall(VecSetBlockSize(Z, bs));
1651: PetscCall(VecPlaceArray(Z, PetscSafePointerPlusOffset(x, start)));
1652: PetscCall(VecRestoreArrayRead(X, &x));
1653: }
1655: /* this is relevant only in debug mode */
1656: PetscCall(VecLockGet(X, &state));
1657: if (state) PetscCall(VecLockReadPush(Z));
1658: Z->ops->placearray = NULL;
1659: Z->ops->replacearray = NULL;
1660: } else { /* Have to create a scatter and do a copy */
1661: PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1662: }
1663: }
1664: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1665: if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1666: PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1667: *Y = Z;
1668: PetscFunctionReturn(PETSC_SUCCESS);
1669: }
1671: /*@
1672: VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`
1674: Collective
1676: Input Parameters:
1677: + X - vector from which subvector was obtained
1678: . is - index set representing the subset of `X`
1679: - Y - subvector being restored
1681: Level: advanced
1683: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1684: @*/
1685: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1686: {
1687: PETSC_UNUSED PetscObjectState dummystate = 0;
1688: PetscBool unchanged;
1690: PetscFunctionBegin;
1693: PetscCheckSameComm(X, 1, is, 2);
1694: PetscAssertPointer(Y, 3);
1697: if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1698: else {
1699: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1700: if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1701: VecScatter scatter;
1702: PetscInt state;
1704: PetscCall(VecLockGet(X, &state));
1705: PetscCheck(state == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec X is locked for read-only or read/write access");
1707: PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1708: if (scatter) {
1709: PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1710: PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1711: } else {
1712: PetscBool iscuda, iship;
1713: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1714: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1716: if (iscuda) {
1717: #if defined(PETSC_HAVE_CUDA)
1718: PetscOffloadMask ymask = (*Y)->offloadmask;
1720: /* The offloadmask of X dictates where to move memory
1721: If X GPU data is valid, then move Y data on GPU if needed
1722: Otherwise, move back to the CPU */
1723: switch (X->offloadmask) {
1724: case PETSC_OFFLOAD_BOTH:
1725: if (ymask == PETSC_OFFLOAD_CPU) {
1726: PetscCall(VecCUDAResetArray(*Y));
1727: } else if (ymask == PETSC_OFFLOAD_GPU) {
1728: X->offloadmask = PETSC_OFFLOAD_GPU;
1729: }
1730: break;
1731: case PETSC_OFFLOAD_GPU:
1732: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1733: break;
1734: case PETSC_OFFLOAD_CPU:
1735: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1736: break;
1737: case PETSC_OFFLOAD_UNALLOCATED:
1738: case PETSC_OFFLOAD_KOKKOS:
1739: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1740: }
1741: #endif
1742: } else if (iship) {
1743: #if defined(PETSC_HAVE_HIP)
1744: PetscOffloadMask ymask = (*Y)->offloadmask;
1746: /* The offloadmask of X dictates where to move memory
1747: If X GPU data is valid, then move Y data on GPU if needed
1748: Otherwise, move back to the CPU */
1749: switch (X->offloadmask) {
1750: case PETSC_OFFLOAD_BOTH:
1751: if (ymask == PETSC_OFFLOAD_CPU) {
1752: PetscCall(VecHIPResetArray(*Y));
1753: } else if (ymask == PETSC_OFFLOAD_GPU) {
1754: X->offloadmask = PETSC_OFFLOAD_GPU;
1755: }
1756: break;
1757: case PETSC_OFFLOAD_GPU:
1758: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1759: break;
1760: case PETSC_OFFLOAD_CPU:
1761: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1762: break;
1763: case PETSC_OFFLOAD_UNALLOCATED:
1764: case PETSC_OFFLOAD_KOKKOS:
1765: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1766: }
1767: #endif
1768: } else {
1769: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1770: PetscCall(VecResetArray(*Y));
1771: }
1772: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1773: }
1774: }
1775: }
1776: PetscCall(VecDestroy(Y));
1777: PetscFunctionReturn(PETSC_SUCCESS);
1778: }
1780: /*@
1781: VecCreateLocalVector - Creates a vector object suitable for use with `VecGetLocalVector()` and friends. You must call `VecDestroy()` when the
1782: vector is no longer needed.
1784: Not Collective.
1786: Input Parameter:
1787: . v - The vector for which the local vector is desired.
1789: Output Parameter:
1790: . w - Upon exit this contains the local vector.
1792: Level: beginner
1794: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1795: @*/
1796: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1797: {
1798: PetscMPIInt size;
1800: PetscFunctionBegin;
1802: PetscAssertPointer(w, 2);
1803: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1804: if (size == 1) PetscCall(VecDuplicate(v, w));
1805: else if (v->ops->createlocalvector) PetscUseTypeMethod(v, createlocalvector, w);
1806: else {
1807: VecType type;
1808: PetscInt n;
1810: PetscCall(VecCreate(PETSC_COMM_SELF, w));
1811: PetscCall(VecGetLocalSize(v, &n));
1812: PetscCall(VecSetSizes(*w, n, n));
1813: PetscCall(VecGetBlockSize(v, &n));
1814: PetscCall(VecSetBlockSize(*w, n));
1815: PetscCall(VecGetType(v, &type));
1816: PetscCall(VecSetType(*w, type));
1817: }
1818: PetscFunctionReturn(PETSC_SUCCESS);
1819: }
1821: /*@
1822: VecGetLocalVectorRead - Maps the local portion of a vector into a
1823: vector.
1825: Not Collective.
1827: Input Parameter:
1828: . v - The vector for which the local vector is desired.
1830: Output Parameter:
1831: . w - Upon exit this contains the local vector.
1833: Level: beginner
1835: Notes:
1836: You must call `VecRestoreLocalVectorRead()` when the local
1837: vector is no longer needed.
1839: This function is similar to `VecGetArrayRead()` which maps the local
1840: portion into a raw pointer. `VecGetLocalVectorRead()` is usually
1841: almost as efficient as `VecGetArrayRead()` but in certain circumstances
1842: `VecGetLocalVectorRead()` can be much more efficient than
1843: `VecGetArrayRead()`. This is because the construction of a contiguous
1844: array representing the vector data required by `VecGetArrayRead()` can
1845: be an expensive operation for certain vector types. For example, for
1846: GPU vectors `VecGetArrayRead()` requires that the data between device
1847: and host is synchronized.
1849: Unlike `VecGetLocalVector()`, this routine is not collective and
1850: preserves cached information.
1852: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1853: @*/
1854: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1855: {
1856: PetscFunctionBegin;
1859: VecCheckSameLocalSize(v, 1, w, 2);
1860: if (v->ops->getlocalvectorread) {
1861: PetscUseTypeMethod(v, getlocalvectorread, w);
1862: } else {
1863: PetscScalar *a;
1865: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1866: PetscCall(VecPlaceArray(w, a));
1867: }
1868: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1869: PetscCall(VecLockReadPush(v));
1870: PetscCall(VecLockReadPush(w));
1871: PetscFunctionReturn(PETSC_SUCCESS);
1872: }
1874: /*@
1875: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1876: previously mapped into a vector using `VecGetLocalVectorRead()`.
1878: Not Collective.
1880: Input Parameters:
1881: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
1882: - w - The vector into which the local portion of `v` was mapped.
1884: Level: beginner
1886: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1887: @*/
1888: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1889: {
1890: PetscFunctionBegin;
1893: if (v->ops->restorelocalvectorread) {
1894: PetscUseTypeMethod(v, restorelocalvectorread, w);
1895: } else {
1896: const PetscScalar *a;
1898: PetscCall(VecGetArrayRead(w, &a));
1899: PetscCall(VecRestoreArrayRead(v, &a));
1900: PetscCall(VecResetArray(w));
1901: }
1902: PetscCall(VecLockReadPop(v));
1903: PetscCall(VecLockReadPop(w));
1904: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1905: PetscFunctionReturn(PETSC_SUCCESS);
1906: }
1908: /*@
1909: VecGetLocalVector - Maps the local portion of a vector into a
1910: vector.
1912: Collective
1914: Input Parameter:
1915: . v - The vector for which the local vector is desired.
1917: Output Parameter:
1918: . w - Upon exit this contains the local vector.
1920: Level: beginner
1922: Notes:
1923: You must call `VecRestoreLocalVector()` when the local
1924: vector is no longer needed.
1926: This function is similar to `VecGetArray()` which maps the local
1927: portion into a raw pointer. `VecGetLocalVector()` is usually about as
1928: efficient as `VecGetArray()` but in certain circumstances
1929: `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1930: This is because the construction of a contiguous array representing
1931: the vector data required by `VecGetArray()` can be an expensive
1932: operation for certain vector types. For example, for GPU vectors
1933: `VecGetArray()` requires that the data between device and host is
1934: synchronized.
1936: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1937: @*/
1938: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1939: {
1940: PetscFunctionBegin;
1943: VecCheckSameLocalSize(v, 1, w, 2);
1944: if (v->ops->getlocalvector) {
1945: PetscUseTypeMethod(v, getlocalvector, w);
1946: } else {
1947: PetscScalar *a;
1949: PetscCall(VecGetArray(v, &a));
1950: PetscCall(VecPlaceArray(w, a));
1951: }
1952: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1953: PetscFunctionReturn(PETSC_SUCCESS);
1954: }
1956: /*@
1957: VecRestoreLocalVector - Unmaps the local portion of a vector
1958: previously mapped into a vector using `VecGetLocalVector()`.
1960: Logically Collective.
1962: Input Parameters:
1963: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
1964: - w - The vector into which the local portion of `v` was mapped.
1966: Level: beginner
1968: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1969: @*/
1970: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1971: {
1972: PetscFunctionBegin;
1975: if (v->ops->restorelocalvector) {
1976: PetscUseTypeMethod(v, restorelocalvector, w);
1977: } else {
1978: PetscScalar *a;
1979: PetscCall(VecGetArray(w, &a));
1980: PetscCall(VecRestoreArray(v, &a));
1981: PetscCall(VecResetArray(w));
1982: }
1983: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1984: PetscCall(PetscObjectStateIncrease((PetscObject)v));
1985: PetscFunctionReturn(PETSC_SUCCESS);
1986: }
1988: /*@C
1989: VecGetArray - Returns a pointer to a contiguous array that contains this
1990: MPI processes's portion of the vector data
1992: Logically Collective
1994: Input Parameter:
1995: . x - the vector
1997: Output Parameter:
1998: . a - location to put pointer to the array
2000: Level: beginner
2002: Notes:
2003: For the standard PETSc vectors, `VecGetArray()` returns a pointer to the local data array and
2004: does not use any copies. If the underlying vector data is not stored in a contiguous array
2005: this routine will copy the data to a contiguous array and return a pointer to that. You MUST
2006: call `VecRestoreArray()` when you no longer need access to the array.
2008: Fortran Notes:
2009: `VecGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayF90()`
2011: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
2012: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2013: @*/
2014: PetscErrorCode VecGetArray(Vec x, PetscScalar **a)
2015: {
2016: PetscFunctionBegin;
2018: PetscCall(VecSetErrorIfLocked(x, 1));
2019: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
2020: PetscUseTypeMethod(x, getarray, a);
2021: } else if (x->petscnative) { /* VECSTANDARD */
2022: *a = *((PetscScalar **)x->data);
2023: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
2024: PetscFunctionReturn(PETSC_SUCCESS);
2025: }
2027: /*@C
2028: VecRestoreArray - Restores a vector after `VecGetArray()` has been called and the array is no longer needed
2030: Logically Collective
2032: Input Parameters:
2033: + x - the vector
2034: - a - location of pointer to array obtained from `VecGetArray()`
2036: Level: beginner
2038: Fortran Notes:
2039: `VecRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayF90()`
2041: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2042: `VecGetArrayPair()`, `VecRestoreArrayPair()`
2043: @*/
2044: PetscErrorCode VecRestoreArray(Vec x, PetscScalar **a)
2045: {
2046: PetscFunctionBegin;
2048: if (a) PetscAssertPointer(a, 2);
2049: if (x->ops->restorearray) {
2050: PetscUseTypeMethod(x, restorearray, a);
2051: } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
2052: if (a) *a = NULL;
2053: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2054: PetscFunctionReturn(PETSC_SUCCESS);
2055: }
2056: /*@C
2057: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
2059: Not Collective
2061: Input Parameter:
2062: . x - the vector
2064: Output Parameter:
2065: . a - the array
2067: Level: beginner
2069: Notes:
2070: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2072: Unlike `VecGetArray()`, preserves cached information like vector norms.
2074: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
2075: implementations may require a copy, but such implementations should cache the contiguous representation so that
2076: only one copy is performed when this routine is called multiple times in sequence.
2078: Fortran Notes:
2079: `VecGetArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayReadF90()`
2081: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2082: @*/
2083: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar **a)
2084: {
2085: PetscFunctionBegin;
2087: PetscAssertPointer(a, 2);
2088: if (x->ops->getarrayread) {
2089: PetscUseTypeMethod(x, getarrayread, a);
2090: } else if (x->ops->getarray) {
2091: PetscObjectState state;
2093: /* VECNEST, VECCUDA, VECKOKKOS etc */
2094: // x->ops->getarray may bump the object state, but since we know this is a read-only get
2095: // we can just undo that
2096: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2097: PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
2098: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2099: } else if (x->petscnative) {
2100: /* VECSTANDARD */
2101: *a = *((PetscScalar **)x->data);
2102: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
2103: PetscFunctionReturn(PETSC_SUCCESS);
2104: }
2106: /*@C
2107: VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`
2109: Not Collective
2111: Input Parameters:
2112: + x - the vector
2113: - a - the array
2115: Level: beginner
2117: Fortran Notes:
2118: `VecRestoreArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayReadF90()`
2120: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2121: @*/
2122: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar **a)
2123: {
2124: PetscFunctionBegin;
2126: if (a) PetscAssertPointer(a, 2);
2127: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
2128: /* nothing */
2129: } else if (x->ops->restorearrayread) { /* VECNEST */
2130: PetscUseTypeMethod(x, restorearrayread, a);
2131: } else { /* No one? */
2132: PetscObjectState state;
2134: // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2135: // we can just undo that
2136: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2137: PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2138: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2139: }
2140: if (a) *a = NULL;
2141: PetscFunctionReturn(PETSC_SUCCESS);
2142: }
2144: /*@C
2145: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
2146: MPI processes's portion of the vector data.
2148: Logically Collective
2150: Input Parameter:
2151: . x - the vector
2153: Output Parameter:
2154: . a - location to put pointer to the array
2156: Level: intermediate
2158: Note:
2159: The values in this array are NOT valid, the caller of this routine is responsible for putting
2160: values into the array; any values it does not set will be invalid.
2162: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2164: For vectors associated with GPUs, the host and device vectors are not synchronized before
2165: giving access. If you need correct values in the array use `VecGetArray()`
2167: Fortran Notes:
2168: `VecGetArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayWriteF90()`
2170: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteF90()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
2171: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
2172: @*/
2173: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar **a)
2174: {
2175: PetscFunctionBegin;
2177: PetscAssertPointer(a, 2);
2178: PetscCall(VecSetErrorIfLocked(x, 1));
2179: if (x->ops->getarraywrite) {
2180: PetscUseTypeMethod(x, getarraywrite, a);
2181: } else {
2182: PetscCall(VecGetArray(x, a));
2183: }
2184: PetscFunctionReturn(PETSC_SUCCESS);
2185: }
2187: /*@C
2188: VecRestoreArrayWrite - Restores a vector after `VecGetArrayWrite()` has been called.
2190: Logically Collective
2192: Input Parameters:
2193: + x - the vector
2194: - a - location of pointer to array obtained from `VecGetArray()`
2196: Level: beginner
2198: Fortran Notes:
2199: `VecRestoreArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayWriteF90()`
2201: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteF90()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2202: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2203: @*/
2204: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar **a)
2205: {
2206: PetscFunctionBegin;
2208: if (a) PetscAssertPointer(a, 2);
2209: if (x->ops->restorearraywrite) {
2210: PetscUseTypeMethod(x, restorearraywrite, a);
2211: } else if (x->ops->restorearray) {
2212: PetscUseTypeMethod(x, restorearray, a);
2213: }
2214: if (a) *a = NULL;
2215: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2216: PetscFunctionReturn(PETSC_SUCCESS);
2217: }
2219: /*@C
2220: VecGetArrays - Returns a pointer to the arrays in a set of vectors
2221: that were created by a call to `VecDuplicateVecs()`.
2223: Logically Collective; No Fortran Support
2225: Input Parameters:
2226: + x - the vectors
2227: - n - the number of vectors
2229: Output Parameter:
2230: . a - location to put pointer to the array
2232: Level: intermediate
2234: Note:
2235: You MUST call `VecRestoreArrays()` when you no longer need access to the arrays.
2237: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2238: @*/
2239: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2240: {
2241: PetscInt i;
2242: PetscScalar **q;
2244: PetscFunctionBegin;
2245: PetscAssertPointer(x, 1);
2247: PetscAssertPointer(a, 3);
2248: PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2249: PetscCall(PetscMalloc1(n, &q));
2250: for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2251: *a = q;
2252: PetscFunctionReturn(PETSC_SUCCESS);
2253: }
2255: /*@C
2256: VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2257: has been called.
2259: Logically Collective; No Fortran Support
2261: Input Parameters:
2262: + x - the vector
2263: . n - the number of vectors
2264: - a - location of pointer to arrays obtained from `VecGetArrays()`
2266: Notes:
2267: For regular PETSc vectors this routine does not involve any copies. For
2268: any special vectors that do not store local vector data in a contiguous
2269: array, this routine will copy the data back into the underlying
2270: vector data structure from the arrays obtained with `VecGetArrays()`.
2272: Level: intermediate
2274: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2275: @*/
2276: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2277: {
2278: PetscInt i;
2279: PetscScalar **q = *a;
2281: PetscFunctionBegin;
2282: PetscAssertPointer(x, 1);
2284: PetscAssertPointer(a, 3);
2286: for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2287: PetscCall(PetscFree(q));
2288: PetscFunctionReturn(PETSC_SUCCESS);
2289: }
2291: /*@C
2292: VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g.,
2293: `VECCUDA`), the returned pointer will be a device pointer to the device memory that contains
2294: this MPI processes's portion of the vector data.
2296: Logically Collective; No Fortran Support
2298: Input Parameter:
2299: . x - the vector
2301: Output Parameters:
2302: + a - location to put pointer to the array
2303: - mtype - memory type of the array
2305: Level: beginner
2307: Note:
2308: Device data is guaranteed to have the latest value. Otherwise, when this is a host vector
2309: (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host
2310: pointer.
2312: For `VECKOKKOS`, if Kokkos is configured without device (e.g., use serial or openmp), per
2313: this function, the vector works like `VECSEQ`/`VECMPI`; otherwise, it works like `VECCUDA` or
2314: `VECHIP` etc.
2316: Use `VecRestoreArrayAndMemType()` when the array access is no longer needed.
2318: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`,
2319: `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2320: @*/
2321: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2322: {
2323: PetscFunctionBegin;
2326: PetscAssertPointer(a, 2);
2327: if (mtype) PetscAssertPointer(mtype, 3);
2328: PetscCall(VecSetErrorIfLocked(x, 1));
2329: if (x->ops->getarrayandmemtype) {
2330: /* VECCUDA, VECKOKKOS etc */
2331: PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2332: } else {
2333: /* VECSTANDARD, VECNEST, VECVIENNACL */
2334: PetscCall(VecGetArray(x, a));
2335: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2336: }
2337: PetscFunctionReturn(PETSC_SUCCESS);
2338: }
2340: /*@C
2341: VecRestoreArrayAndMemType - Restores a vector after `VecGetArrayAndMemType()` has been called.
2343: Logically Collective; No Fortran Support
2345: Input Parameters:
2346: + x - the vector
2347: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`
2349: Level: beginner
2351: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`,
2352: `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2353: @*/
2354: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar **a)
2355: {
2356: PetscFunctionBegin;
2359: if (a) PetscAssertPointer(a, 2);
2360: if (x->ops->restorearrayandmemtype) {
2361: /* VECCUDA, VECKOKKOS etc */
2362: PetscUseTypeMethod(x, restorearrayandmemtype, a);
2363: } else {
2364: /* VECNEST, VECVIENNACL */
2365: PetscCall(VecRestoreArray(x, a));
2366: } /* VECSTANDARD does nothing */
2367: if (a) *a = NULL;
2368: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2369: PetscFunctionReturn(PETSC_SUCCESS);
2370: }
2372: /*@C
2373: VecGetArrayReadAndMemType - Like `VecGetArrayRead()`, but if the input vector is a device vector, it will return a read-only device pointer.
2374: The returned pointer is guaranteed to point to up-to-date data. For host vectors, it functions as `VecGetArrayRead()`.
2376: Not Collective; No Fortran Support
2378: Input Parameter:
2379: . x - the vector
2381: Output Parameters:
2382: + a - the array
2383: - mtype - memory type of the array
2385: Level: beginner
2387: Notes:
2388: The array must be returned using a matching call to `VecRestoreArrayReadAndMemType()`.
2390: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2391: @*/
2392: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar **a, PetscMemType *mtype)
2393: {
2394: PetscFunctionBegin;
2397: PetscAssertPointer(a, 2);
2398: if (mtype) PetscAssertPointer(mtype, 3);
2399: if (x->ops->getarrayreadandmemtype) {
2400: /* VECCUDA/VECHIP though they are also petscnative */
2401: PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2402: } else if (x->ops->getarrayandmemtype) {
2403: /* VECKOKKOS */
2404: PetscObjectState state;
2406: // see VecGetArrayRead() for why
2407: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2408: PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2409: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2410: } else {
2411: PetscCall(VecGetArrayRead(x, a));
2412: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2413: }
2414: PetscFunctionReturn(PETSC_SUCCESS);
2415: }
2417: /*@C
2418: VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`
2420: Not Collective; No Fortran Support
2422: Input Parameters:
2423: + x - the vector
2424: - a - the array
2426: Level: beginner
2428: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2429: @*/
2430: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar **a)
2431: {
2432: PetscFunctionBegin;
2435: if (a) PetscAssertPointer(a, 2);
2436: if (x->ops->restorearrayreadandmemtype) {
2437: /* VECCUDA/VECHIP */
2438: PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2439: } else if (!x->petscnative) {
2440: /* VECNEST */
2441: PetscCall(VecRestoreArrayRead(x, a));
2442: }
2443: if (a) *a = NULL;
2444: PetscFunctionReturn(PETSC_SUCCESS);
2445: }
2447: /*@C
2448: VecGetArrayWriteAndMemType - Like `VecGetArrayWrite()`, but if this is a device vector it will always return
2449: a device pointer to the device memory that contains this processor's portion of the vector data.
2451: Logically Collective; No Fortran Support
2453: Input Parameter:
2454: . x - the vector
2456: Output Parameters:
2457: + a - the array
2458: - mtype - memory type of the array
2460: Level: beginner
2462: Note:
2463: The array must be returned using a matching call to `VecRestoreArrayWriteAndMemType()`, where it will label the device memory as most recent.
2465: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2466: @*/
2467: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2468: {
2469: PetscFunctionBegin;
2472: PetscCall(VecSetErrorIfLocked(x, 1));
2473: PetscAssertPointer(a, 2);
2474: if (mtype) PetscAssertPointer(mtype, 3);
2475: if (x->ops->getarraywriteandmemtype) {
2476: /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2477: PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2478: } else if (x->ops->getarrayandmemtype) {
2479: PetscCall(VecGetArrayAndMemType(x, a, mtype));
2480: } else {
2481: /* VECNEST, VECVIENNACL */
2482: PetscCall(VecGetArrayWrite(x, a));
2483: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2484: }
2485: PetscFunctionReturn(PETSC_SUCCESS);
2486: }
2488: /*@C
2489: VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`
2491: Logically Collective; No Fortran Support
2493: Input Parameters:
2494: + x - the vector
2495: - a - the array
2497: Level: beginner
2499: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2500: @*/
2501: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar **a)
2502: {
2503: PetscFunctionBegin;
2506: PetscCall(VecSetErrorIfLocked(x, 1));
2507: if (a) PetscAssertPointer(a, 2);
2508: if (x->ops->restorearraywriteandmemtype) {
2509: /* VECCUDA/VECHIP */
2510: PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2511: PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2512: } else if (x->ops->restorearrayandmemtype) {
2513: PetscCall(VecRestoreArrayAndMemType(x, a));
2514: } else {
2515: PetscCall(VecRestoreArray(x, a));
2516: }
2517: if (a) *a = NULL;
2518: PetscFunctionReturn(PETSC_SUCCESS);
2519: }
2521: /*@
2522: VecPlaceArray - Allows one to replace the array in a vector with an
2523: array provided by the user. This is useful to avoid copying an array
2524: into a vector.
2526: Logically Collective; No Fortran Support
2528: Input Parameters:
2529: + vec - the vector
2530: - array - the array
2532: Level: developer
2534: Notes:
2535: Use `VecReplaceArray()` instead to permanently replace the array
2537: You can return to the original array with a call to `VecResetArray()`. `vec` does not take
2538: ownership of `array` in any way.
2540: The user must free `array` themselves but be careful not to
2541: do so before the vector has either been destroyed, had its original array restored with
2542: `VecResetArray()` or permanently replaced with `VecReplaceArray()`.
2544: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2545: @*/
2546: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2547: {
2548: PetscFunctionBegin;
2551: if (array) PetscAssertPointer(array, 2);
2552: PetscUseTypeMethod(vec, placearray, array);
2553: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2554: PetscFunctionReturn(PETSC_SUCCESS);
2555: }
2557: /*@C
2558: VecReplaceArray - Allows one to replace the array in a vector with an
2559: array provided by the user. This is useful to avoid copying an array
2560: into a vector.
2562: Logically Collective; No Fortran Support
2564: Input Parameters:
2565: + vec - the vector
2566: - array - the array
2568: Level: developer
2570: Notes:
2571: This permanently replaces the array and frees the memory associated
2572: with the old array. Use `VecPlaceArray()` to temporarily replace the array.
2574: The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
2575: freed by the user. It will be freed when the vector is destroyed.
2577: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2578: @*/
2579: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2580: {
2581: PetscFunctionBegin;
2584: PetscUseTypeMethod(vec, replacearray, array);
2585: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2586: PetscFunctionReturn(PETSC_SUCCESS);
2587: }
2589: /*MC
2590: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
2591: and makes them accessible via a Fortran pointer.
2593: Synopsis:
2594: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
2596: Collective
2598: Input Parameters:
2599: + x - a vector to mimic
2600: - n - the number of vectors to obtain
2602: Output Parameters:
2603: + y - Fortran pointer to the array of vectors
2604: - ierr - error code
2606: Example of Usage:
2607: .vb
2608: #include <petsc/finclude/petscvec.h>
2609: use petscvec
2611: Vec x
2612: Vec, pointer :: y(:)
2613: ....
2614: call VecDuplicateVecsF90(x,2,y,ierr)
2615: call VecSet(y(2),alpha,ierr)
2616: call VecSet(y(2),alpha,ierr)
2617: ....
2618: call VecDestroyVecsF90(2,y,ierr)
2619: .ve
2621: Level: beginner
2623: Note:
2624: Use `VecDestroyVecsF90()` to free the space.
2626: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecsF90()`, `VecDuplicateVecs()`
2627: M*/
2629: /*MC
2630: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2631: `VecGetArrayF90()`.
2633: Synopsis:
2634: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2636: Logically Collective
2638: Input Parameters:
2639: + x - vector
2640: - xx_v - the Fortran pointer to the array
2642: Output Parameter:
2643: . ierr - error code
2645: Example of Usage:
2646: .vb
2647: #include <petsc/finclude/petscvec.h>
2648: use petscvec
2650: PetscScalar, pointer :: xx_v(:)
2651: ....
2652: call VecGetArrayF90(x,xx_v,ierr)
2653: xx_v(3) = a
2654: call VecRestoreArrayF90(x,xx_v,ierr)
2655: .ve
2657: Level: beginner
2659: .seealso: [](ch_vectors), `Vec`, `VecGetArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrayReadF90()`
2660: M*/
2662: /*MC
2663: VecDestroyVecsF90 - Frees a block of vectors obtained with `VecDuplicateVecsF90()`.
2665: Synopsis:
2666: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
2668: Collective
2670: Input Parameters:
2671: + n - the number of vectors previously obtained
2672: - x - pointer to array of vector pointers
2674: Output Parameter:
2675: . ierr - error code
2677: Level: beginner
2679: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecs()`, `VecDuplicateVecsF90()`
2680: M*/
2682: /*MC
2683: VecGetArrayF90 - Accesses a vector array from Fortran. For default PETSc
2684: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2685: this routine is implementation dependent. You MUST call `VecRestoreArrayF90()`
2686: when you no longer need access to the array.
2688: Synopsis:
2689: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2691: Logically Collective
2693: Input Parameter:
2694: . x - vector
2696: Output Parameters:
2697: + xx_v - the Fortran pointer to the array
2698: - ierr - error code
2700: Example of Usage:
2701: .vb
2702: #include <petsc/finclude/petscvec.h>
2703: use petscvec
2705: PetscScalar, pointer :: xx_v(:)
2706: ....
2707: call VecGetArrayF90(x,xx_v,ierr)
2708: xx_v(3) = a
2709: call VecRestoreArrayF90(x,xx_v,ierr)
2710: .ve
2712: Level: beginner
2714: Note:
2715: If you ONLY intend to read entries from the array and not change any entries you should use `VecGetArrayReadF90()`.
2717: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayReadF90()`
2718: M*/
2720: /*MC
2721: VecGetArrayReadF90 - Accesses a read only array from Fortran. For default PETSc
2722: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2723: this routine is implementation dependent. You MUST call `VecRestoreArrayReadF90()`
2724: when you no longer need access to the array.
2726: Synopsis:
2727: VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2729: Logically Collective
2731: Input Parameter:
2732: . x - vector
2734: Output Parameters:
2735: + xx_v - the Fortran pointer to the array
2736: - ierr - error code
2738: Example of Usage:
2739: .vb
2740: #include <petsc/finclude/petscvec.h>
2741: use petscvec
2743: PetscScalar, pointer :: xx_v(:)
2744: ....
2745: call VecGetArrayReadF90(x,xx_v,ierr)
2746: a = xx_v(3)
2747: call VecRestoreArrayReadF90(x,xx_v,ierr)
2748: .ve
2750: Level: beginner
2752: Note:
2753: If you intend to write entries into the array you must use `VecGetArrayF90()`.
2755: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecGetArrayF90()`
2756: M*/
2758: /*MC
2759: VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2760: `VecGetArrayReadF90()`.
2762: Synopsis:
2763: VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2765: Logically Collective
2767: Input Parameters:
2768: + x - vector
2769: - xx_v - the Fortran pointer to the array
2771: Output Parameter:
2772: . ierr - error code
2774: Example of Usage:
2775: .vb
2776: #include <petsc/finclude/petscvec.h>
2777: use petscvec
2779: PetscScalar, pointer :: xx_v(:)
2780: ....
2781: call VecGetArrayReadF90(x,xx_v,ierr)
2782: a = xx_v(3)
2783: call VecRestoreArrayReadF90(x,xx_v,ierr)
2784: .ve
2786: Level: beginner
2788: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecRestoreArrayF90()`
2789: M*/
2791: /*@C
2792: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2793: processor's portion of the vector data. You MUST call `VecRestoreArray2d()`
2794: when you no longer need access to the array.
2796: Logically Collective
2798: Input Parameters:
2799: + x - the vector
2800: . m - first dimension of two dimensional array
2801: . n - second dimension of two dimensional array
2802: . mstart - first index you will use in first coordinate direction (often 0)
2803: - nstart - first index in the second coordinate direction (often 0)
2805: Output Parameter:
2806: . a - location to put pointer to the array
2808: Level: developer
2810: Notes:
2811: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2812: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2813: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2814: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2816: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2818: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2819: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2820: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2821: @*/
2822: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2823: {
2824: PetscInt i, N;
2825: PetscScalar *aa;
2827: PetscFunctionBegin;
2829: PetscAssertPointer(a, 6);
2831: PetscCall(VecGetLocalSize(x, &N));
2832: 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);
2833: PetscCall(VecGetArray(x, &aa));
2835: PetscCall(PetscMalloc1(m, a));
2836: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2837: *a -= mstart;
2838: PetscFunctionReturn(PETSC_SUCCESS);
2839: }
2841: /*@C
2842: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2843: processor's portion of the vector data. You MUST call `VecRestoreArray2dWrite()`
2844: when you no longer need access to the array.
2846: Logically Collective
2848: Input Parameters:
2849: + x - the vector
2850: . m - first dimension of two dimensional array
2851: . n - second dimension of two dimensional array
2852: . mstart - first index you will use in first coordinate direction (often 0)
2853: - nstart - first index in the second coordinate direction (often 0)
2855: Output Parameter:
2856: . a - location to put pointer to the array
2858: Level: developer
2860: Notes:
2861: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2862: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2863: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2864: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2866: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2868: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2869: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2870: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2871: @*/
2872: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2873: {
2874: PetscInt i, N;
2875: PetscScalar *aa;
2877: PetscFunctionBegin;
2879: PetscAssertPointer(a, 6);
2881: PetscCall(VecGetLocalSize(x, &N));
2882: 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);
2883: PetscCall(VecGetArrayWrite(x, &aa));
2885: PetscCall(PetscMalloc1(m, a));
2886: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2887: *a -= mstart;
2888: PetscFunctionReturn(PETSC_SUCCESS);
2889: }
2891: /*@C
2892: VecRestoreArray2d - Restores a vector after `VecGetArray2d()` has been called.
2894: Logically Collective
2896: Input Parameters:
2897: + x - the vector
2898: . m - first dimension of two dimensional array
2899: . n - second dimension of the two dimensional array
2900: . mstart - first index you will use in first coordinate direction (often 0)
2901: . nstart - first index in the second coordinate direction (often 0)
2902: - a - location of pointer to array obtained from `VecGetArray2d()`
2904: Level: developer
2906: Notes:
2907: For regular PETSc vectors this routine does not involve any copies. For
2908: any special vectors that do not store local vector data in a contiguous
2909: array, this routine will copy the data back into the underlying
2910: vector data structure from the array obtained with `VecGetArray()`.
2912: This routine actually zeros out the `a` pointer.
2914: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2915: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2916: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2917: @*/
2918: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2919: {
2920: void *dummy;
2922: PetscFunctionBegin;
2924: PetscAssertPointer(a, 6);
2926: dummy = (void *)(*a + mstart);
2927: PetscCall(PetscFree(dummy));
2928: PetscCall(VecRestoreArray(x, NULL));
2929: *a = NULL;
2930: PetscFunctionReturn(PETSC_SUCCESS);
2931: }
2933: /*@C
2934: VecRestoreArray2dWrite - Restores a vector after `VecGetArray2dWrite()` has been called.
2936: Logically Collective
2938: Input Parameters:
2939: + x - the vector
2940: . m - first dimension of two dimensional array
2941: . n - second dimension of the two dimensional array
2942: . mstart - first index you will use in first coordinate direction (often 0)
2943: . nstart - first index in the second coordinate direction (often 0)
2944: - a - location of pointer to array obtained from `VecGetArray2d()`
2946: Level: developer
2948: Notes:
2949: For regular PETSc vectors this routine does not involve any copies. For
2950: any special vectors that do not store local vector data in a contiguous
2951: array, this routine will copy the data back into the underlying
2952: vector data structure from the array obtained with `VecGetArray()`.
2954: This routine actually zeros out the `a` pointer.
2956: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2957: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2958: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2959: @*/
2960: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2961: {
2962: void *dummy;
2964: PetscFunctionBegin;
2966: PetscAssertPointer(a, 6);
2968: dummy = (void *)(*a + mstart);
2969: PetscCall(PetscFree(dummy));
2970: PetscCall(VecRestoreArrayWrite(x, NULL));
2971: PetscFunctionReturn(PETSC_SUCCESS);
2972: }
2974: /*@C
2975: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2976: processor's portion of the vector data. You MUST call `VecRestoreArray1d()`
2977: when you no longer need access to the array.
2979: Logically Collective
2981: Input Parameters:
2982: + x - the vector
2983: . m - first dimension of two dimensional array
2984: - mstart - first index you will use in first coordinate direction (often 0)
2986: Output Parameter:
2987: . a - location to put pointer to the array
2989: Level: developer
2991: Notes:
2992: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2993: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2994: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2996: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2998: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2999: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3000: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3001: @*/
3002: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3003: {
3004: PetscInt N;
3006: PetscFunctionBegin;
3008: PetscAssertPointer(a, 4);
3010: PetscCall(VecGetLocalSize(x, &N));
3011: 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);
3012: PetscCall(VecGetArray(x, a));
3013: *a -= mstart;
3014: PetscFunctionReturn(PETSC_SUCCESS);
3015: }
3017: /*@C
3018: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
3019: processor's portion of the vector data. You MUST call `VecRestoreArray1dWrite()`
3020: when you no longer need access to the array.
3022: Logically Collective
3024: Input Parameters:
3025: + x - the vector
3026: . m - first dimension of two dimensional array
3027: - mstart - first index you will use in first coordinate direction (often 0)
3029: Output Parameter:
3030: . a - location to put pointer to the array
3032: Level: developer
3034: Notes:
3035: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3036: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3037: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3039: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3041: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3042: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3043: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3044: @*/
3045: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3046: {
3047: PetscInt N;
3049: PetscFunctionBegin;
3051: PetscAssertPointer(a, 4);
3053: PetscCall(VecGetLocalSize(x, &N));
3054: 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);
3055: PetscCall(VecGetArrayWrite(x, a));
3056: *a -= mstart;
3057: PetscFunctionReturn(PETSC_SUCCESS);
3058: }
3060: /*@C
3061: VecRestoreArray1d - Restores a vector after `VecGetArray1d()` has been called.
3063: Logically Collective
3065: Input Parameters:
3066: + x - the vector
3067: . m - first dimension of two dimensional array
3068: . mstart - first index you will use in first coordinate direction (often 0)
3069: - a - location of pointer to array obtained from `VecGetArray1d()`
3071: Level: developer
3073: Notes:
3074: For regular PETSc vectors this routine does not involve any copies. For
3075: any special vectors that do not store local vector data in a contiguous
3076: array, this routine will copy the data back into the underlying
3077: vector data structure from the array obtained with `VecGetArray1d()`.
3079: This routine actually zeros out the `a` pointer.
3081: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3082: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3083: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3084: @*/
3085: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3086: {
3087: PetscFunctionBegin;
3090: PetscCall(VecRestoreArray(x, NULL));
3091: *a = NULL;
3092: PetscFunctionReturn(PETSC_SUCCESS);
3093: }
3095: /*@C
3096: VecRestoreArray1dWrite - Restores a vector after `VecGetArray1dWrite()` has been called.
3098: Logically Collective
3100: Input Parameters:
3101: + x - the vector
3102: . m - first dimension of two dimensional array
3103: . mstart - first index you will use in first coordinate direction (often 0)
3104: - a - location of pointer to array obtained from `VecGetArray1d()`
3106: Level: developer
3108: Notes:
3109: For regular PETSc vectors this routine does not involve any copies. For
3110: any special vectors that do not store local vector data in a contiguous
3111: array, this routine will copy the data back into the underlying
3112: vector data structure from the array obtained with `VecGetArray1d()`.
3114: This routine actually zeros out the `a` pointer.
3116: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3117: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3118: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3119: @*/
3120: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3121: {
3122: PetscFunctionBegin;
3125: PetscCall(VecRestoreArrayWrite(x, NULL));
3126: *a = NULL;
3127: PetscFunctionReturn(PETSC_SUCCESS);
3128: }
3130: /*@C
3131: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
3132: processor's portion of the vector data. You MUST call `VecRestoreArray3d()`
3133: when you no longer need access to the array.
3135: Logically Collective
3137: Input Parameters:
3138: + x - the vector
3139: . m - first dimension of three dimensional array
3140: . n - second dimension of three dimensional array
3141: . p - third dimension of three dimensional array
3142: . mstart - first index you will use in first coordinate direction (often 0)
3143: . nstart - first index in the second coordinate direction (often 0)
3144: - pstart - first index in the third coordinate direction (often 0)
3146: Output Parameter:
3147: . a - location to put pointer to the array
3149: Level: developer
3151: Notes:
3152: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3153: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3154: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3155: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3157: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3159: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3160: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
3161: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3162: @*/
3163: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3164: {
3165: PetscInt i, N, j;
3166: PetscScalar *aa, **b;
3168: PetscFunctionBegin;
3170: PetscAssertPointer(a, 8);
3172: PetscCall(VecGetLocalSize(x, &N));
3173: 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);
3174: PetscCall(VecGetArray(x, &aa));
3176: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3177: b = (PetscScalar **)((*a) + m);
3178: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3179: for (i = 0; i < m; i++)
3180: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset(aa, i * n * p + j * p - pstart);
3181: *a -= mstart;
3182: PetscFunctionReturn(PETSC_SUCCESS);
3183: }
3185: /*@C
3186: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
3187: processor's portion of the vector data. You MUST call `VecRestoreArray3dWrite()`
3188: when you no longer need access to the array.
3190: Logically Collective
3192: Input Parameters:
3193: + x - the vector
3194: . m - first dimension of three dimensional array
3195: . n - second dimension of three dimensional array
3196: . p - third dimension of three dimensional array
3197: . mstart - first index you will use in first coordinate direction (often 0)
3198: . nstart - first index in the second coordinate direction (often 0)
3199: - pstart - first index in the third coordinate direction (often 0)
3201: Output Parameter:
3202: . a - location to put pointer to the array
3204: Level: developer
3206: Notes:
3207: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3208: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3209: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3210: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3212: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3214: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3215: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3216: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3217: @*/
3218: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3219: {
3220: PetscInt i, N, j;
3221: PetscScalar *aa, **b;
3223: PetscFunctionBegin;
3225: PetscAssertPointer(a, 8);
3227: PetscCall(VecGetLocalSize(x, &N));
3228: 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);
3229: PetscCall(VecGetArrayWrite(x, &aa));
3231: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3232: b = (PetscScalar **)((*a) + m);
3233: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3234: for (i = 0; i < m; i++)
3235: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3237: *a -= mstart;
3238: PetscFunctionReturn(PETSC_SUCCESS);
3239: }
3241: /*@C
3242: VecRestoreArray3d - Restores a vector after `VecGetArray3d()` has been called.
3244: Logically Collective
3246: Input Parameters:
3247: + x - the vector
3248: . m - first dimension of three dimensional array
3249: . n - second dimension of the three dimensional array
3250: . p - third dimension of the three dimensional array
3251: . mstart - first index you will use in first coordinate direction (often 0)
3252: . nstart - first index in the second coordinate direction (often 0)
3253: . pstart - first index in the third coordinate direction (often 0)
3254: - a - location of pointer to array obtained from VecGetArray3d()
3256: Level: developer
3258: Notes:
3259: For regular PETSc vectors this routine does not involve any copies. For
3260: any special vectors that do not store local vector data in a contiguous
3261: array, this routine will copy the data back into the underlying
3262: vector data structure from the array obtained with `VecGetArray()`.
3264: This routine actually zeros out the `a` pointer.
3266: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3267: `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3268: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3269: @*/
3270: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3271: {
3272: void *dummy;
3274: PetscFunctionBegin;
3276: PetscAssertPointer(a, 8);
3278: dummy = (void *)(*a + mstart);
3279: PetscCall(PetscFree(dummy));
3280: PetscCall(VecRestoreArray(x, NULL));
3281: *a = NULL;
3282: PetscFunctionReturn(PETSC_SUCCESS);
3283: }
3285: /*@C
3286: VecRestoreArray3dWrite - Restores a vector after `VecGetArray3dWrite()` has been called.
3288: Logically Collective
3290: Input Parameters:
3291: + x - the vector
3292: . m - first dimension of three dimensional array
3293: . n - second dimension of the three dimensional array
3294: . p - third dimension of the three dimensional array
3295: . mstart - first index you will use in first coordinate direction (often 0)
3296: . nstart - first index in the second coordinate direction (often 0)
3297: . pstart - first index in the third coordinate direction (often 0)
3298: - a - location of pointer to array obtained from VecGetArray3d()
3300: Level: developer
3302: Notes:
3303: For regular PETSc vectors this routine does not involve any copies. For
3304: any special vectors that do not store local vector data in a contiguous
3305: array, this routine will copy the data back into the underlying
3306: vector data structure from the array obtained with `VecGetArray()`.
3308: This routine actually zeros out the `a` pointer.
3310: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3311: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3312: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3313: @*/
3314: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3315: {
3316: void *dummy;
3318: PetscFunctionBegin;
3320: PetscAssertPointer(a, 8);
3322: dummy = (void *)(*a + mstart);
3323: PetscCall(PetscFree(dummy));
3324: PetscCall(VecRestoreArrayWrite(x, NULL));
3325: *a = NULL;
3326: PetscFunctionReturn(PETSC_SUCCESS);
3327: }
3329: /*@C
3330: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3331: processor's portion of the vector data. You MUST call `VecRestoreArray4d()`
3332: when you no longer need access to the array.
3334: Logically Collective
3336: Input Parameters:
3337: + x - the vector
3338: . m - first dimension of four dimensional array
3339: . n - second dimension of four dimensional array
3340: . p - third dimension of four dimensional array
3341: . q - fourth dimension of four dimensional array
3342: . mstart - first index you will use in first coordinate direction (often 0)
3343: . nstart - first index in the second coordinate direction (often 0)
3344: . pstart - first index in the third coordinate direction (often 0)
3345: - qstart - first index in the fourth coordinate direction (often 0)
3347: Output Parameter:
3348: . a - location to put pointer to the array
3350: Level: developer
3352: Notes:
3353: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3354: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3355: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3356: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3358: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3360: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3361: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3362: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3363: @*/
3364: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3365: {
3366: PetscInt i, N, j, k;
3367: PetscScalar *aa, ***b, **c;
3369: PetscFunctionBegin;
3371: PetscAssertPointer(a, 10);
3373: PetscCall(VecGetLocalSize(x, &N));
3374: 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);
3375: PetscCall(VecGetArray(x, &aa));
3377: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3378: b = (PetscScalar ***)((*a) + m);
3379: c = (PetscScalar **)(b + m * n);
3380: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3381: for (i = 0; i < m; i++)
3382: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3383: for (i = 0; i < m; i++)
3384: for (j = 0; j < n; j++)
3385: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3386: *a -= mstart;
3387: PetscFunctionReturn(PETSC_SUCCESS);
3388: }
3390: /*@C
3391: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3392: processor's portion of the vector data. You MUST call `VecRestoreArray4dWrite()`
3393: when you no longer need access to the array.
3395: Logically Collective
3397: Input Parameters:
3398: + x - the vector
3399: . m - first dimension of four dimensional array
3400: . n - second dimension of four dimensional array
3401: . p - third dimension of four dimensional array
3402: . q - fourth dimension of four dimensional array
3403: . mstart - first index you will use in first coordinate direction (often 0)
3404: . nstart - first index in the second coordinate direction (often 0)
3405: . pstart - first index in the third coordinate direction (often 0)
3406: - qstart - first index in the fourth coordinate direction (often 0)
3408: Output Parameter:
3409: . a - location to put pointer to the array
3411: Level: developer
3413: Notes:
3414: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3415: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3416: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3417: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3419: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3421: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3422: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3423: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3424: @*/
3425: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3426: {
3427: PetscInt i, N, j, k;
3428: PetscScalar *aa, ***b, **c;
3430: PetscFunctionBegin;
3432: PetscAssertPointer(a, 10);
3434: PetscCall(VecGetLocalSize(x, &N));
3435: 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);
3436: PetscCall(VecGetArrayWrite(x, &aa));
3438: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3439: b = (PetscScalar ***)((*a) + m);
3440: c = (PetscScalar **)(b + m * n);
3441: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3442: for (i = 0; i < m; i++)
3443: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3444: for (i = 0; i < m; i++)
3445: for (j = 0; j < n; j++)
3446: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3447: *a -= mstart;
3448: PetscFunctionReturn(PETSC_SUCCESS);
3449: }
3451: /*@C
3452: VecRestoreArray4d - Restores a vector after `VecGetArray4d()` has been called.
3454: Logically Collective
3456: Input Parameters:
3457: + x - the vector
3458: . m - first dimension of four dimensional array
3459: . n - second dimension of the four dimensional array
3460: . p - third dimension of the four dimensional array
3461: . q - fourth dimension of the four dimensional array
3462: . mstart - first index you will use in first coordinate direction (often 0)
3463: . nstart - first index in the second coordinate direction (often 0)
3464: . pstart - first index in the third coordinate direction (often 0)
3465: . qstart - first index in the fourth coordinate direction (often 0)
3466: - a - location of pointer to array obtained from VecGetArray4d()
3468: Level: developer
3470: Notes:
3471: For regular PETSc vectors this routine does not involve any copies. For
3472: any special vectors that do not store local vector data in a contiguous
3473: array, this routine will copy the data back into the underlying
3474: vector data structure from the array obtained with `VecGetArray()`.
3476: This routine actually zeros out the `a` pointer.
3478: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3479: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3480: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecGet`
3481: @*/
3482: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3483: {
3484: void *dummy;
3486: PetscFunctionBegin;
3488: PetscAssertPointer(a, 10);
3490: dummy = (void *)(*a + mstart);
3491: PetscCall(PetscFree(dummy));
3492: PetscCall(VecRestoreArray(x, NULL));
3493: *a = NULL;
3494: PetscFunctionReturn(PETSC_SUCCESS);
3495: }
3497: /*@C
3498: VecRestoreArray4dWrite - Restores a vector after `VecGetArray4dWrite()` has been called.
3500: Logically Collective
3502: Input Parameters:
3503: + x - the vector
3504: . m - first dimension of four dimensional array
3505: . n - second dimension of the four dimensional array
3506: . p - third dimension of the four dimensional array
3507: . q - fourth dimension of the four dimensional array
3508: . mstart - first index you will use in first coordinate direction (often 0)
3509: . nstart - first index in the second coordinate direction (often 0)
3510: . pstart - first index in the third coordinate direction (often 0)
3511: . qstart - first index in the fourth coordinate direction (often 0)
3512: - a - location of pointer to array obtained from `VecGetArray4d()`
3514: Level: developer
3516: Notes:
3517: For regular PETSc vectors this routine does not involve any copies. For
3518: any special vectors that do not store local vector data in a contiguous
3519: array, this routine will copy the data back into the underlying
3520: vector data structure from the array obtained with `VecGetArray()`.
3522: This routine actually zeros out the `a` pointer.
3524: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3525: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3526: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3527: @*/
3528: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3529: {
3530: void *dummy;
3532: PetscFunctionBegin;
3534: PetscAssertPointer(a, 10);
3536: dummy = (void *)(*a + mstart);
3537: PetscCall(PetscFree(dummy));
3538: PetscCall(VecRestoreArrayWrite(x, NULL));
3539: *a = NULL;
3540: PetscFunctionReturn(PETSC_SUCCESS);
3541: }
3543: /*@C
3544: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3545: processor's portion of the vector data. You MUST call `VecRestoreArray2dRead()`
3546: when you no longer need access to the array.
3548: Logically Collective
3550: Input Parameters:
3551: + x - the vector
3552: . m - first dimension of two dimensional array
3553: . n - second dimension of two dimensional array
3554: . mstart - first index you will use in first coordinate direction (often 0)
3555: - nstart - first index in the second coordinate direction (often 0)
3557: Output Parameter:
3558: . a - location to put pointer to the array
3560: Level: developer
3562: Notes:
3563: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
3564: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3565: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3566: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
3568: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3570: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3571: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3572: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3573: @*/
3574: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3575: {
3576: PetscInt i, N;
3577: const PetscScalar *aa;
3579: PetscFunctionBegin;
3581: PetscAssertPointer(a, 6);
3583: PetscCall(VecGetLocalSize(x, &N));
3584: 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);
3585: PetscCall(VecGetArrayRead(x, &aa));
3587: PetscCall(PetscMalloc1(m, a));
3588: for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3589: *a -= mstart;
3590: PetscFunctionReturn(PETSC_SUCCESS);
3591: }
3593: /*@C
3594: VecRestoreArray2dRead - Restores a vector after `VecGetArray2dRead()` has been called.
3596: Logically Collective
3598: Input Parameters:
3599: + x - the vector
3600: . m - first dimension of two dimensional array
3601: . n - second dimension of the two dimensional array
3602: . mstart - first index you will use in first coordinate direction (often 0)
3603: . nstart - first index in the second coordinate direction (often 0)
3604: - a - location of pointer to array obtained from VecGetArray2d()
3606: Level: developer
3608: Notes:
3609: For regular PETSc vectors this routine does not involve any copies. For
3610: any special vectors that do not store local vector data in a contiguous
3611: array, this routine will copy the data back into the underlying
3612: vector data structure from the array obtained with `VecGetArray()`.
3614: This routine actually zeros out the `a` pointer.
3616: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3617: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3618: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3619: @*/
3620: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3621: {
3622: void *dummy;
3624: PetscFunctionBegin;
3626: PetscAssertPointer(a, 6);
3628: dummy = (void *)(*a + mstart);
3629: PetscCall(PetscFree(dummy));
3630: PetscCall(VecRestoreArrayRead(x, NULL));
3631: *a = NULL;
3632: PetscFunctionReturn(PETSC_SUCCESS);
3633: }
3635: /*@C
3636: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3637: processor's portion of the vector data. You MUST call `VecRestoreArray1dRead()`
3638: when you no longer need access to the array.
3640: Logically Collective
3642: Input Parameters:
3643: + x - the vector
3644: . m - first dimension of two dimensional array
3645: - mstart - first index you will use in first coordinate direction (often 0)
3647: Output Parameter:
3648: . a - location to put pointer to the array
3650: Level: developer
3652: Notes:
3653: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3654: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3655: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3657: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3659: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3660: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3661: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3662: @*/
3663: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3664: {
3665: PetscInt N;
3667: PetscFunctionBegin;
3669: PetscAssertPointer(a, 4);
3671: PetscCall(VecGetLocalSize(x, &N));
3672: 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);
3673: PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3674: *a -= mstart;
3675: PetscFunctionReturn(PETSC_SUCCESS);
3676: }
3678: /*@C
3679: VecRestoreArray1dRead - Restores a vector after `VecGetArray1dRead()` has been called.
3681: Logically Collective
3683: Input Parameters:
3684: + x - the vector
3685: . m - first dimension of two dimensional array
3686: . mstart - first index you will use in first coordinate direction (often 0)
3687: - a - location of pointer to array obtained from `VecGetArray1dRead()`
3689: Level: developer
3691: Notes:
3692: For regular PETSc vectors this routine does not involve any copies. For
3693: any special vectors that do not store local vector data in a contiguous
3694: array, this routine will copy the data back into the underlying
3695: vector data structure from the array obtained with `VecGetArray1dRead()`.
3697: This routine actually zeros out the `a` pointer.
3699: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3700: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3701: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3702: @*/
3703: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3704: {
3705: PetscFunctionBegin;
3708: PetscCall(VecRestoreArrayRead(x, NULL));
3709: *a = NULL;
3710: PetscFunctionReturn(PETSC_SUCCESS);
3711: }
3713: /*@C
3714: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3715: processor's portion of the vector data. You MUST call `VecRestoreArray3dRead()`
3716: when you no longer need access to the array.
3718: Logically Collective
3720: Input Parameters:
3721: + x - the vector
3722: . m - first dimension of three dimensional array
3723: . n - second dimension of three dimensional array
3724: . p - third dimension of three dimensional array
3725: . mstart - first index you will use in first coordinate direction (often 0)
3726: . nstart - first index in the second coordinate direction (often 0)
3727: - pstart - first index in the third coordinate direction (often 0)
3729: Output Parameter:
3730: . a - location to put pointer to the array
3732: Level: developer
3734: Notes:
3735: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3736: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3737: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3738: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3dRead()`.
3740: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3742: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3743: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3744: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3745: @*/
3746: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3747: {
3748: PetscInt i, N, j;
3749: const PetscScalar *aa;
3750: PetscScalar **b;
3752: PetscFunctionBegin;
3754: PetscAssertPointer(a, 8);
3756: PetscCall(VecGetLocalSize(x, &N));
3757: 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);
3758: PetscCall(VecGetArrayRead(x, &aa));
3760: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3761: b = (PetscScalar **)((*a) + m);
3762: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3763: for (i = 0; i < m; i++)
3764: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset((PetscScalar *)aa, i * n * p + j * p - pstart);
3765: *a -= mstart;
3766: PetscFunctionReturn(PETSC_SUCCESS);
3767: }
3769: /*@C
3770: VecRestoreArray3dRead - Restores a vector after `VecGetArray3dRead()` has been called.
3772: Logically Collective
3774: Input Parameters:
3775: + x - the vector
3776: . m - first dimension of three dimensional array
3777: . n - second dimension of the three dimensional array
3778: . p - third dimension of the three dimensional array
3779: . mstart - first index you will use in first coordinate direction (often 0)
3780: . nstart - first index in the second coordinate direction (often 0)
3781: . pstart - first index in the third coordinate direction (often 0)
3782: - a - location of pointer to array obtained from `VecGetArray3dRead()`
3784: Level: developer
3786: Notes:
3787: For regular PETSc vectors this routine does not involve any copies. For
3788: any special vectors that do not store local vector data in a contiguous
3789: array, this routine will copy the data back into the underlying
3790: vector data structure from the array obtained with `VecGetArray()`.
3792: This routine actually zeros out the `a` pointer.
3794: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3795: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3796: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3797: @*/
3798: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3799: {
3800: void *dummy;
3802: PetscFunctionBegin;
3804: PetscAssertPointer(a, 8);
3806: dummy = (void *)(*a + mstart);
3807: PetscCall(PetscFree(dummy));
3808: PetscCall(VecRestoreArrayRead(x, NULL));
3809: *a = NULL;
3810: PetscFunctionReturn(PETSC_SUCCESS);
3811: }
3813: /*@C
3814: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3815: processor's portion of the vector data. You MUST call `VecRestoreArray4dRead()`
3816: when you no longer need access to the array.
3818: Logically Collective
3820: Input Parameters:
3821: + x - the vector
3822: . m - first dimension of four dimensional array
3823: . n - second dimension of four dimensional array
3824: . p - third dimension of four dimensional array
3825: . q - fourth dimension of four dimensional array
3826: . mstart - first index you will use in first coordinate direction (often 0)
3827: . nstart - first index in the second coordinate direction (often 0)
3828: . pstart - first index in the third coordinate direction (often 0)
3829: - qstart - first index in the fourth coordinate direction (often 0)
3831: Output Parameter:
3832: . a - location to put pointer to the array
3834: Level: beginner
3836: Notes:
3837: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3838: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3839: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3840: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3842: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3844: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3845: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3846: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3847: @*/
3848: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3849: {
3850: PetscInt i, N, j, k;
3851: const PetscScalar *aa;
3852: PetscScalar ***b, **c;
3854: PetscFunctionBegin;
3856: PetscAssertPointer(a, 10);
3858: PetscCall(VecGetLocalSize(x, &N));
3859: 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);
3860: PetscCall(VecGetArrayRead(x, &aa));
3862: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3863: b = (PetscScalar ***)((*a) + m);
3864: c = (PetscScalar **)(b + m * n);
3865: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3866: for (i = 0; i < m; i++)
3867: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3868: for (i = 0; i < m; i++)
3869: for (j = 0; j < n; j++)
3870: 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;
3871: *a -= mstart;
3872: PetscFunctionReturn(PETSC_SUCCESS);
3873: }
3875: /*@C
3876: VecRestoreArray4dRead - Restores a vector after `VecGetArray4d()` has been called.
3878: Logically Collective
3880: Input Parameters:
3881: + x - the vector
3882: . m - first dimension of four dimensional array
3883: . n - second dimension of the four dimensional array
3884: . p - third dimension of the four dimensional array
3885: . q - fourth dimension of the four dimensional array
3886: . mstart - first index you will use in first coordinate direction (often 0)
3887: . nstart - first index in the second coordinate direction (often 0)
3888: . pstart - first index in the third coordinate direction (often 0)
3889: . qstart - first index in the fourth coordinate direction (often 0)
3890: - a - location of pointer to array obtained from `VecGetArray4dRead()`
3892: Level: beginner
3894: Notes:
3895: For regular PETSc vectors this routine does not involve any copies. For
3896: any special vectors that do not store local vector data in a contiguous
3897: array, this routine will copy the data back into the underlying
3898: vector data structure from the array obtained with `VecGetArray()`.
3900: This routine actually zeros out the `a` pointer.
3902: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3903: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3904: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3905: @*/
3906: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3907: {
3908: void *dummy;
3910: PetscFunctionBegin;
3912: PetscAssertPointer(a, 10);
3914: dummy = (void *)(*a + mstart);
3915: PetscCall(PetscFree(dummy));
3916: PetscCall(VecRestoreArrayRead(x, NULL));
3917: *a = NULL;
3918: PetscFunctionReturn(PETSC_SUCCESS);
3919: }
3921: /*@
3922: VecLockGet - Get the current lock status of a vector
3924: Logically Collective
3926: Input Parameter:
3927: . x - the vector
3929: Output Parameter:
3930: . state - greater than zero indicates the vector is locked for read; less than zero indicates the vector is
3931: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3933: Level: advanced
3935: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3936: @*/
3937: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3938: {
3939: PetscFunctionBegin;
3941: PetscAssertPointer(state, 2);
3942: *state = x->lock;
3943: PetscFunctionReturn(PETSC_SUCCESS);
3944: }
3946: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3947: {
3948: PetscFunctionBegin;
3950: PetscAssertPointer(file, 2);
3951: PetscAssertPointer(func, 3);
3952: PetscAssertPointer(line, 4);
3953: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3954: {
3955: const int index = x->lockstack.currentsize - 1;
3957: *file = index < 0 ? NULL : x->lockstack.file[index];
3958: *func = index < 0 ? NULL : x->lockstack.function[index];
3959: *line = index < 0 ? 0 : x->lockstack.line[index];
3960: }
3961: #else
3962: *file = NULL;
3963: *func = NULL;
3964: *line = 0;
3965: #endif
3966: PetscFunctionReturn(PETSC_SUCCESS);
3967: }
3969: /*@
3970: VecLockReadPush - Push a read-only lock on a vector to prevent it from being written to
3972: Logically Collective
3974: Input Parameter:
3975: . x - the vector
3977: Level: intermediate
3979: Notes:
3980: If this is set then calls to `VecGetArray()` or `VecSetValues()` or any other routines that change the vectors values will generate an error.
3982: The call can be nested, i.e., called multiple times on the same vector, but each `VecLockReadPush()` has to have one matching
3983: `VecLockReadPop()`, which removes the latest read-only lock.
3985: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3986: @*/
3987: PetscErrorCode VecLockReadPush(Vec x)
3988: {
3989: PetscFunctionBegin;
3991: 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");
3992: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3993: {
3994: const char *file, *func;
3995: int index, line;
3997: if ((index = petscstack.currentsize - 2) < 0) {
3998: // vec was locked "outside" of petsc, either in user-land or main. the error message will
3999: // now show this function as the culprit, but it will include the stacktrace
4000: file = "unknown user-file";
4001: func = "unknown_user_function";
4002: line = 0;
4003: } else {
4004: file = petscstack.file[index];
4005: func = petscstack.function[index];
4006: line = petscstack.line[index];
4007: }
4008: PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
4009: }
4010: #endif
4011: PetscFunctionReturn(PETSC_SUCCESS);
4012: }
4014: /*@
4015: VecLockReadPop - Pop a read-only lock from a vector
4017: Logically Collective
4019: Input Parameter:
4020: . x - the vector
4022: Level: intermediate
4024: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
4025: @*/
4026: PetscErrorCode VecLockReadPop(Vec x)
4027: {
4028: PetscFunctionBegin;
4030: PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
4031: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
4032: {
4033: const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];
4035: PetscStackPop_Private(x->lockstack, previous);
4036: }
4037: #endif
4038: PetscFunctionReturn(PETSC_SUCCESS);
4039: }
4041: /*@
4042: VecLockWriteSet - Lock or unlock a vector for exclusive read/write access
4044: Logically Collective
4046: Input Parameters:
4047: + x - the vector
4048: - flg - `PETSC_TRUE` to lock the vector for exclusive read/write access; `PETSC_FALSE` to unlock it.
4050: Level: intermediate
4052: Notes:
4053: The function is useful in split-phase computations, which usually have a begin phase and an end phase.
4054: One can call `VecLockWriteSet`(x,`PETSC_TRUE`) in the begin phase to lock a vector for exclusive
4055: access, and call `VecLockWriteSet`(x,`PETSC_FALSE`) in the end phase to unlock the vector from exclusive
4056: access. In this way, one is ensured no other operations can access the vector in between. The code may like
4058: .vb
4059: VecGetArray(x,&xdata); // begin phase
4060: VecLockWriteSet(v,PETSC_TRUE);
4062: Other operations, which can not access x anymore (they can access xdata, of course)
4064: VecRestoreArray(x,&vdata); // end phase
4065: VecLockWriteSet(v,PETSC_FALSE);
4066: .ve
4068: The call can not be nested on the same vector, in other words, one can not call `VecLockWriteSet`(x,`PETSC_TRUE`)
4069: again before calling `VecLockWriteSet`(v,`PETSC_FALSE`).
4071: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
4072: @*/
4073: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
4074: {
4075: PetscFunctionBegin;
4077: if (flg) {
4078: 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");
4079: 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");
4080: x->lock = -1;
4081: } else {
4082: 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");
4083: x->lock = 0;
4084: }
4085: PetscFunctionReturn(PETSC_SUCCESS);
4086: }