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: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
332: @*/
333: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
334: {
335: PetscFunctionBegin;
338: VecCheckAssembled(x);
339: if (p) PetscAssertPointer(p, 2);
340: PetscAssertPointer(val, 3);
341: PetscCall(VecLockReadPush(x));
342: PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
343: PetscUseTypeMethod(x, max, p, val);
344: PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
345: PetscCall(VecLockReadPop(x));
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: /*@C
350: VecMin - Determines the vector component with minimum real part and its location.
352: Collective
354: Input Parameter:
355: . x - the vector
357: Output Parameters:
358: + p - the index of `val` (pass `NULL` if you don't want this location) in the vector
359: - val - the minimum component
361: Level: intermediate
363: Notes:
364: Returns the value `PETSC_MAX_REAL` and negative `p` if the vector is of length 0.
366: This returns the smallest index with the minimum value
368: .seealso: [](ch_vectors), `Vec`, `VecMax()`
369: @*/
370: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
371: {
372: PetscFunctionBegin;
375: VecCheckAssembled(x);
376: if (p) PetscAssertPointer(p, 2);
377: PetscAssertPointer(val, 3);
378: PetscCall(VecLockReadPush(x));
379: PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
380: PetscUseTypeMethod(x, min, p, val);
381: PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
382: PetscCall(VecLockReadPop(x));
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: /*@
387: VecTDot - Computes an indefinite vector dot product. That is, this
388: routine does NOT use the complex conjugate.
390: Collective
392: Input Parameters:
393: + x - first vector
394: - y - second vector
396: Output Parameter:
397: . val - the dot product
399: Level: intermediate
401: Notes for Users of Complex Numbers:
402: For complex vectors, `VecTDot()` computes the indefinite form
403: $ val = (x,y) = y^T x,
404: where y^T denotes the transpose of y.
406: Use `VecDot()` for the inner product
407: $ val = (x,y) = y^H x,
408: where y^H denotes the conjugate transpose of y.
410: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
411: @*/
412: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
413: {
414: PetscFunctionBegin;
417: PetscAssertPointer(val, 3);
420: PetscCheckSameTypeAndComm(x, 1, y, 2);
421: VecCheckSameSize(x, 1, y, 2);
422: VecCheckAssembled(x);
423: VecCheckAssembled(y);
425: PetscCall(VecLockReadPush(x));
426: PetscCall(VecLockReadPush(y));
427: PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
428: PetscUseTypeMethod(x, tdot, y, val);
429: PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
430: PetscCall(VecLockReadPop(x));
431: PetscCall(VecLockReadPop(y));
432: PetscFunctionReturn(PETSC_SUCCESS);
433: }
435: PetscErrorCode VecScaleAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
436: {
437: PetscReal norms[4];
438: PetscBool flgs[4];
439: PetscScalar one = 1.0;
441: PetscFunctionBegin;
444: VecCheckAssembled(x);
445: PetscCall(VecSetErrorIfLocked(x, 1));
447: if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);
449: /* get current stashed norms */
450: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
452: PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
453: VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
454: PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));
456: PetscCall(PetscObjectStateIncrease((PetscObject)x));
457: /* put the scaled stashed norms back into the Vec */
458: for (PetscInt i = 0; i < 4; i++) {
459: PetscReal ar = PetscAbsScalar(alpha);
460: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
461: }
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }
465: /*@
466: VecScale - Scales a vector.
468: Logically Collective
470: Input Parameters:
471: + x - the vector
472: - alpha - the scalar
474: Level: intermediate
476: Note:
477: For a vector with n components, `VecScale()` computes x[i] = alpha * x[i], for i=1,...,n.
479: .seealso: [](ch_vectors), `Vec`, `VecSet()`
480: @*/
481: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
482: {
483: PetscFunctionBegin;
484: PetscCall(VecScaleAsync_Private(x, alpha, NULL));
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
489: {
490: PetscFunctionBegin;
493: VecCheckAssembled(x);
495: PetscCall(VecSetErrorIfLocked(x, 1));
497: if (alpha == 0) {
498: PetscReal norm;
499: PetscBool set;
501: PetscCall(VecNormAvailable(x, NORM_2, &set, &norm));
502: if (set == PETSC_TRUE && norm == 0) PetscFunctionReturn(PETSC_SUCCESS);
503: }
504: PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
505: VecMethodDispatch(x, dctx, VecAsyncFnName(Set), set, (Vec, PetscScalar, PetscDeviceContext), alpha);
506: PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
507: PetscCall(PetscObjectStateIncrease((PetscObject)x));
509: /* norms can be simply set (if |alpha|*N not too large) */
510: {
511: PetscReal val = PetscAbsScalar(alpha);
512: const PetscInt N = x->map->N;
514: if (N == 0) {
515: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0));
516: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
517: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
518: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
519: } else if (val > PETSC_MAX_REAL / N) {
520: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
521: } else {
522: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
523: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
524: val *= PetscSqrtReal((PetscReal)N);
525: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
526: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
527: }
528: }
529: PetscFunctionReturn(PETSC_SUCCESS);
530: }
532: /*@
533: VecSet - Sets all components of a vector to a single scalar value.
535: Logically Collective
537: Input Parameters:
538: + x - the vector
539: - alpha - the scalar
541: Level: beginner
543: Notes:
544: For a vector of dimension n, `VecSet()` sets x[i] = alpha, for i=1,...,n,
545: so that all vector entries then equal the identical
546: scalar value, `alpha`. Use the more general routine
547: `VecSetValues()` to set different vector entries.
549: You CANNOT call this after you have called `VecSetValues()` but before you call
550: `VecAssemblyBegin()`
552: If `alpha` is zero and the norm of the vector is known to be zero then this skips the unneeded zeroing process
554: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
555: @*/
556: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
557: {
558: PetscFunctionBegin;
559: PetscCall(VecSetAsync_Private(x, alpha, NULL));
560: PetscFunctionReturn(PETSC_SUCCESS);
561: }
563: PetscErrorCode VecAXPYAsync_Private(Vec y, PetscScalar alpha, Vec x, PetscDeviceContext dctx)
564: {
565: PetscFunctionBegin;
570: PetscCheckSameTypeAndComm(x, 3, y, 1);
571: VecCheckSameSize(x, 3, y, 1);
572: VecCheckAssembled(x);
573: VecCheckAssembled(y);
575: if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
576: PetscCall(VecSetErrorIfLocked(y, 1));
577: if (x == y) {
578: PetscCall(VecScale(y, alpha + 1.0));
579: PetscFunctionReturn(PETSC_SUCCESS);
580: }
581: PetscCall(VecLockReadPush(x));
582: PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
583: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPY), axpy, (Vec, PetscScalar, Vec, PetscDeviceContext), alpha, x);
584: PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
585: PetscCall(VecLockReadPop(x));
586: PetscCall(PetscObjectStateIncrease((PetscObject)y));
587: PetscFunctionReturn(PETSC_SUCCESS);
588: }
589: /*@
590: VecAXPY - Computes `y = alpha x + y`.
592: Logically Collective
594: Input Parameters:
595: + alpha - the scalar
596: . x - vector scale by `alpha`
597: - y - vector accumulated into
599: Output Parameter:
600: . y - output vector
602: Level: intermediate
604: Notes:
605: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
606: .vb
607: VecAXPY(y,alpha,x) y = alpha x + y
608: VecAYPX(y,beta,x) y = x + beta y
609: VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
610: VecWAXPY(w,alpha,x,y) w = alpha x + y
611: VecAXPBYPCZ(z,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
612: VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
613: .ve
615: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
616: @*/
617: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
618: {
619: PetscFunctionBegin;
620: PetscCall(VecAXPYAsync_Private(y, alpha, x, NULL));
621: PetscFunctionReturn(PETSC_SUCCESS);
622: }
624: PetscErrorCode VecAYPXAsync_Private(Vec y, PetscScalar beta, Vec x, PetscDeviceContext dctx)
625: {
626: PetscFunctionBegin;
631: PetscCheckSameTypeAndComm(x, 3, y, 1);
632: VecCheckSameSize(x, 1, y, 3);
633: VecCheckAssembled(x);
634: VecCheckAssembled(y);
636: PetscCall(VecSetErrorIfLocked(y, 1));
637: if (x == y) {
638: PetscCall(VecScale(y, beta + 1.0));
639: PetscFunctionReturn(PETSC_SUCCESS);
640: }
641: PetscCall(VecLockReadPush(x));
642: if (beta == (PetscScalar)0.0) {
643: PetscCall(VecCopy(x, y));
644: } else {
645: PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
646: VecMethodDispatch(y, dctx, VecAsyncFnName(AYPX), aypx, (Vec, PetscScalar, Vec, PetscDeviceContext), beta, x);
647: PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
648: PetscCall(PetscObjectStateIncrease((PetscObject)y));
649: }
650: PetscCall(VecLockReadPop(x));
651: PetscFunctionReturn(PETSC_SUCCESS);
652: }
654: /*@
655: VecAYPX - Computes `y = x + beta y`.
657: Logically Collective
659: Input Parameters:
660: + beta - the scalar
661: . x - the unscaled vector
662: - y - the vector to be scaled
664: Output Parameter:
665: . y - output vector
667: Level: intermediate
669: Developer Notes:
670: The implementation is optimized for `beta` of -1.0, 0.0, and 1.0
672: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
673: @*/
674: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
675: {
676: PetscFunctionBegin;
677: PetscCall(VecAYPXAsync_Private(y, beta, x, NULL));
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: PetscErrorCode VecAXPBYAsync_Private(Vec y, PetscScalar alpha, PetscScalar beta, Vec x, PetscDeviceContext dctx)
682: {
683: PetscFunctionBegin;
688: PetscCheckSameTypeAndComm(x, 4, y, 1);
689: VecCheckSameSize(y, 1, x, 4);
690: VecCheckAssembled(x);
691: VecCheckAssembled(y);
694: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
695: if (x == y) {
696: PetscCall(VecScale(y, alpha + beta));
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: PetscCall(VecSetErrorIfLocked(y, 1));
701: PetscCall(VecLockReadPush(x));
702: PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
703: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPBY), axpby, (Vec, PetscScalar, PetscScalar, Vec, PetscDeviceContext), alpha, beta, x);
704: PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
705: PetscCall(PetscObjectStateIncrease((PetscObject)y));
706: PetscCall(VecLockReadPop(x));
707: PetscFunctionReturn(PETSC_SUCCESS);
708: }
710: /*@
711: VecAXPBY - Computes `y = alpha x + beta y`.
713: Logically Collective
715: Input Parameters:
716: + alpha - first scalar
717: . beta - second scalar
718: . x - the first scaled vector
719: - y - the second scaled vector
721: Output Parameter:
722: . y - output vector
724: Level: intermediate
726: Developer Notes:
727: The implementation is optimized for `alpha` and/or `beta` values of 0.0 and 1.0
729: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
730: @*/
731: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
732: {
733: PetscFunctionBegin;
734: PetscCall(VecAXPBYAsync_Private(y, alpha, beta, x, NULL));
735: PetscFunctionReturn(PETSC_SUCCESS);
736: }
738: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
739: {
740: PetscFunctionBegin;
747: PetscCheckSameTypeAndComm(x, 5, y, 6);
748: PetscCheckSameTypeAndComm(x, 5, z, 1);
749: VecCheckSameSize(x, 5, y, 6);
750: VecCheckSameSize(x, 5, z, 1);
751: PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
752: PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
753: VecCheckAssembled(x);
754: VecCheckAssembled(y);
755: VecCheckAssembled(z);
759: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
761: PetscCall(VecSetErrorIfLocked(z, 1));
762: PetscCall(VecLockReadPush(x));
763: PetscCall(VecLockReadPush(y));
764: PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
765: VecMethodDispatch(z, dctx, VecAsyncFnName(AXPBYPCZ), axpbypcz, (Vec, PetscScalar, PetscScalar, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, beta, gamma, x, y);
766: PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
767: PetscCall(PetscObjectStateIncrease((PetscObject)z));
768: PetscCall(VecLockReadPop(x));
769: PetscCall(VecLockReadPop(y));
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
772: /*@
773: VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`
775: Logically Collective
777: Input Parameters:
778: + alpha - first scalar
779: . beta - second scalar
780: . gamma - third scalar
781: . x - first vector
782: . y - second vector
783: - z - third vector
785: Output Parameter:
786: . z - output vector
788: Level: intermediate
790: Note:
791: `x`, `y` and `z` must be different vectors
793: Developer Notes:
794: The implementation is optimized for `alpha` of 1.0 and `gamma` of 1.0 or 0.0
796: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
797: @*/
798: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
799: {
800: PetscFunctionBegin;
801: PetscCall(VecAXPBYPCZAsync_Private(z, alpha, beta, gamma, x, y, NULL));
802: PetscFunctionReturn(PETSC_SUCCESS);
803: }
805: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
806: {
807: PetscFunctionBegin;
814: PetscCheckSameTypeAndComm(x, 3, y, 4);
815: PetscCheckSameTypeAndComm(y, 4, w, 1);
816: VecCheckSameSize(x, 3, y, 4);
817: VecCheckSameSize(x, 3, w, 1);
818: PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
819: PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
820: VecCheckAssembled(x);
821: VecCheckAssembled(y);
823: PetscCall(VecSetErrorIfLocked(w, 1));
825: PetscCall(VecLockReadPush(x));
826: PetscCall(VecLockReadPush(y));
827: if (alpha == (PetscScalar)0.0) {
828: PetscCall(VecCopyAsync_Private(y, w, dctx));
829: } else {
830: PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
831: VecMethodDispatch(w, dctx, VecAsyncFnName(WAXPY), waxpy, (Vec, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, x, y);
832: PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
833: PetscCall(PetscObjectStateIncrease((PetscObject)w));
834: }
835: PetscCall(VecLockReadPop(x));
836: PetscCall(VecLockReadPop(y));
837: PetscFunctionReturn(PETSC_SUCCESS);
838: }
840: /*@
841: VecWAXPY - Computes `w = alpha x + y`.
843: Logically Collective
845: Input Parameters:
846: + alpha - the scalar
847: . x - first vector, multiplied by `alpha`
848: - y - second vector
850: Output Parameter:
851: . w - the result
853: Level: intermediate
855: Note:
856: `w` cannot be either `x` or `y`, but `x` and `y` can be the same
858: Developer Notes:
859: The implementation is optimized for alpha of -1.0, 0.0, and 1.0
861: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
862: @*/
863: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
864: {
865: PetscFunctionBegin;
866: PetscCall(VecWAXPYAsync_Private(w, alpha, x, y, NULL));
867: PetscFunctionReturn(PETSC_SUCCESS);
868: }
870: /*@C
871: VecSetValues - Inserts or adds values into certain locations of a vector.
873: Not Collective
875: Input Parameters:
876: + x - vector to insert in
877: . ni - number of elements to add
878: . ix - indices where to add
879: . y - array of values
880: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries
882: Level: beginner
884: Notes:
885: .vb
886: `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
887: .ve
889: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
890: options cannot be mixed without intervening calls to the assembly
891: routines.
893: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
894: MUST be called after all calls to `VecSetValues()` have been completed.
896: VecSetValues() uses 0-based indices in Fortran as well as in C.
898: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
899: negative indices may be passed in ix. These rows are
900: simply ignored. This allows easily inserting element load matrices
901: with homogeneous Dirichlet boundary conditions that you don't want represented
902: in the vector.
904: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
905: `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
906: @*/
907: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
908: {
909: PetscFunctionBeginHot;
911: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
912: PetscAssertPointer(ix, 3);
913: PetscAssertPointer(y, 4);
916: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
917: PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
918: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
919: PetscCall(PetscObjectStateIncrease((PetscObject)x));
920: PetscFunctionReturn(PETSC_SUCCESS);
921: }
923: /*@C
924: VecGetValues - Gets values from certain locations of a vector. Currently
925: can only get values on the same processor on which they are owned
927: Not Collective
929: Input Parameters:
930: + x - vector to get values from
931: . ni - number of elements to get
932: - ix - indices where to get them from (in global 1d numbering)
934: Output Parameter:
935: . y - array of values
937: Level: beginner
939: Notes:
940: The user provides the allocated array y; it is NOT allocated in this routine
942: `VecGetValues()` gets y[i] = x[ix[i]], for i=0,...,ni-1.
944: `VecAssemblyBegin()` and `VecAssemblyEnd()` MUST be called before calling this if `VecSetValues()` or related routine has been called
946: VecGetValues() uses 0-based indices in Fortran as well as in C.
948: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
949: negative indices may be passed in ix. These rows are
950: simply ignored.
952: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
953: @*/
954: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
955: {
956: PetscFunctionBegin;
958: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
959: PetscAssertPointer(ix, 3);
960: PetscAssertPointer(y, 4);
962: VecCheckAssembled(x);
963: PetscUseTypeMethod(x, getvalues, ni, ix, y);
964: PetscFunctionReturn(PETSC_SUCCESS);
965: }
967: /*@C
968: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
970: Not Collective
972: Input Parameters:
973: + x - vector to insert in
974: . ni - number of blocks to add
975: . ix - indices where to add in block count, rather than element count
976: . y - array of values
977: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries
979: Level: intermediate
981: Notes:
982: `VecSetValuesBlocked()` sets x[bs*ix[i]+j] = y[bs*i+j],
983: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
985: Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
986: options cannot be mixed without intervening calls to the assembly
987: routines.
989: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
990: MUST be called after all calls to `VecSetValuesBlocked()` have been completed.
992: `VecSetValuesBlocked()` uses 0-based indices in Fortran as well as in C.
994: Negative indices may be passed in ix, these rows are
995: simply ignored. This allows easily inserting element load matrices
996: with homogeneous Dirichlet boundary conditions that you don't want represented
997: in the vector.
999: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
1000: `VecSetValues()`
1001: @*/
1002: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1003: {
1004: PetscFunctionBeginHot;
1006: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1007: PetscAssertPointer(ix, 3);
1008: PetscAssertPointer(y, 4);
1011: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1012: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1013: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1014: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1015: PetscFunctionReturn(PETSC_SUCCESS);
1016: }
1018: /*@C
1019: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1020: using a local ordering of the nodes.
1022: Not Collective
1024: Input Parameters:
1025: + x - vector to insert in
1026: . ni - number of elements to add
1027: . ix - indices where to add
1028: . y - array of values
1029: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1031: Level: intermediate
1033: Notes:
1034: `VecSetValuesLocal()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
1036: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
1037: options cannot be mixed without intervening calls to the assembly
1038: routines.
1040: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1041: MUST be called after all calls to `VecSetValuesLocal()` have been completed.
1043: `VecSetValuesLocal()` uses 0-based indices in Fortran as well as in C.
1045: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1046: `VecSetValuesBlockedLocal()`
1047: @*/
1048: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1049: {
1050: PetscInt lixp[128], *lix = lixp;
1052: PetscFunctionBeginHot;
1054: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1055: PetscAssertPointer(ix, 3);
1056: PetscAssertPointer(y, 4);
1059: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1060: if (!x->ops->setvalueslocal) {
1061: if (x->map->mapping) {
1062: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1063: PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1064: PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1065: if (ni > 128) PetscCall(PetscFree(lix));
1066: } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1067: } else PetscUseTypeMethod(x, setvalueslocal, ni, ix, y, iora);
1068: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1069: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1070: PetscFunctionReturn(PETSC_SUCCESS);
1071: }
1073: /*@
1074: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1075: using a local ordering of the nodes.
1077: Not Collective
1079: Input Parameters:
1080: + x - vector to insert in
1081: . ni - number of blocks to add
1082: . ix - indices where to add in block count, not element count
1083: . y - array of values
1084: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1086: Level: intermediate
1088: Notes:
1089: `VecSetValuesBlockedLocal()` sets x[bs*ix[i]+j] = y[bs*i+j],
1090: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with `VecSetBlockSize()`.
1092: Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1093: options cannot be mixed without intervening calls to the assembly
1094: routines.
1096: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1097: MUST be called after all calls to `VecSetValuesBlockedLocal()` have been completed.
1099: `VecSetValuesBlockedLocal()` uses 0-based indices in Fortran as well as in C.
1101: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1102: `VecSetLocalToGlobalMapping()`
1103: @*/
1104: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1105: {
1106: PetscInt lixp[128], *lix = lixp;
1108: PetscFunctionBeginHot;
1110: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1111: PetscAssertPointer(ix, 3);
1112: PetscAssertPointer(y, 4);
1114: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1115: if (x->map->mapping) {
1116: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1117: PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1118: PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1119: if (ni > 128) PetscCall(PetscFree(lix));
1120: } else {
1121: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1122: }
1123: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1124: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1125: PetscFunctionReturn(PETSC_SUCCESS);
1126: }
1128: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1129: {
1130: PetscFunctionBegin;
1133: VecCheckAssembled(x);
1135: if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1136: PetscAssertPointer(y, 3);
1137: for (PetscInt i = 0; i < nv; ++i) {
1140: PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1141: VecCheckSameSize(x, 1, y[i], 3);
1142: VecCheckAssembled(y[i]);
1143: PetscCall(VecLockReadPush(y[i]));
1144: }
1145: PetscAssertPointer(result, 4);
1148: PetscCall(VecLockReadPush(x));
1149: PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1150: PetscCall((*mxdot)(x, nv, y, result));
1151: PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1152: PetscCall(VecLockReadPop(x));
1153: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1154: PetscFunctionReturn(PETSC_SUCCESS);
1155: }
1157: /*@
1158: VecMTDot - Computes indefinite vector multiple dot products.
1159: That is, it does NOT use the complex conjugate.
1161: Collective
1163: Input Parameters:
1164: + x - one vector
1165: . nv - number of vectors
1166: - y - array of vectors. Note that vectors are pointers
1168: Output Parameter:
1169: . val - array of the dot products
1171: Level: intermediate
1173: Notes for Users of Complex Numbers:
1174: For complex vectors, `VecMTDot()` computes the indefinite form
1175: $ val = (x,y) = y^T x,
1176: where y^T denotes the transpose of y.
1178: Use `VecMDot()` for the inner product
1179: $ val = (x,y) = y^H x,
1180: where y^H denotes the conjugate transpose of y.
1182: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1183: @*/
1184: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1185: {
1186: PetscFunctionBegin;
1188: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1189: PetscFunctionReturn(PETSC_SUCCESS);
1190: }
1192: /*@
1193: VecMDot - Computes multiple vector dot products.
1195: Collective
1197: Input Parameters:
1198: + x - one vector
1199: . nv - number of vectors
1200: - y - array of vectors.
1202: Output Parameter:
1203: . val - array of the dot products (does not allocate the array)
1205: Level: intermediate
1207: Notes for Users of Complex Numbers:
1208: For complex vectors, `VecMDot()` computes
1209: $ val = (x,y) = y^H x,
1210: where y^H denotes the conjugate transpose of y.
1212: Use `VecMTDot()` for the indefinite form
1213: $ val = (x,y) = y^T x,
1214: where y^T denotes the transpose of y.
1216: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`
1217: @*/
1218: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1219: {
1220: PetscFunctionBegin;
1222: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1223: PetscFunctionReturn(PETSC_SUCCESS);
1224: }
1226: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
1227: {
1228: PetscFunctionBegin;
1230: VecCheckAssembled(y);
1232: PetscCall(VecSetErrorIfLocked(y, 1));
1233: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1234: if (nv) {
1235: PetscInt zeros = 0;
1237: PetscAssertPointer(alpha, 3);
1238: PetscAssertPointer(x, 4);
1239: for (PetscInt i = 0; i < nv; ++i) {
1243: PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1244: VecCheckSameSize(y, 1, x[i], 4);
1245: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1246: VecCheckAssembled(x[i]);
1247: PetscCall(VecLockReadPush(x[i]));
1248: zeros += alpha[i] == (PetscScalar)0.0;
1249: }
1251: if (zeros < nv) {
1252: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1253: VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1254: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1255: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1256: }
1258: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1259: }
1260: PetscFunctionReturn(PETSC_SUCCESS);
1261: }
1263: /*@
1264: VecMAXPY - Computes `y = y + sum alpha[i] x[i]`
1266: Logically Collective
1268: Input Parameters:
1269: + nv - number of scalars and x-vectors
1270: . alpha - array of scalars
1271: . y - one vector
1272: - x - array of vectors
1274: Level: intermediate
1276: Note:
1277: `y` cannot be any of the `x` vectors
1279: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`,`VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1280: @*/
1281: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1282: {
1283: PetscFunctionBegin;
1284: PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1285: PetscFunctionReturn(PETSC_SUCCESS);
1286: }
1288: /*@
1289: VecMAXPBY - Computes `y = beta y + sum alpha[i] x[i]`
1291: Logically Collective
1293: Input Parameters:
1294: + nv - number of scalars and x-vectors
1295: . alpha - array of scalars
1296: . beta - scalar
1297: . y - one vector
1298: - x - array of vectors
1300: Level: intermediate
1302: Note:
1303: `y` cannot be any of the `x` vectors.
1305: Developer Notes:
1306: This is a convenience routine, but implementations might be able to optimize it, for example, when `beta` is zero.
1308: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1309: @*/
1310: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
1311: {
1312: PetscFunctionBegin;
1314: VecCheckAssembled(y);
1316: PetscCall(VecSetErrorIfLocked(y, 1));
1317: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1320: if (y->ops->maxpby) {
1321: PetscInt zeros = 0;
1323: if (nv) {
1324: PetscAssertPointer(alpha, 3);
1325: PetscAssertPointer(x, 5);
1326: }
1328: for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
1332: PetscCheckSameTypeAndComm(y, 1, x[i], 5);
1333: VecCheckSameSize(y, 1, x[i], 5);
1334: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1335: VecCheckAssembled(x[i]);
1336: PetscCall(VecLockReadPush(x[i]));
1337: zeros += alpha[i] == (PetscScalar)0.0;
1338: }
1340: if (zeros < nv) { // has nonzero alpha
1341: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1342: PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1343: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1344: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1345: } else {
1346: PetscCall(VecScale(y, beta));
1347: }
1349: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1350: } else { // no maxpby
1351: if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1352: else PetscCall(VecScale(y, beta));
1353: PetscCall(VecMAXPY(y, nv, alpha, x));
1354: }
1355: PetscFunctionReturn(PETSC_SUCCESS);
1356: }
1358: /*@
1359: VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
1360: in the order they appear in the array. The concatenated vector resides on the same
1361: communicator and is the same type as the source vectors.
1363: Collective
1365: Input Parameters:
1366: + nx - number of vectors to be concatenated
1367: - X - array containing the vectors to be concatenated in the order of concatenation
1369: Output Parameters:
1370: + Y - concatenated vector
1371: - x_is - array of index sets corresponding to the concatenated components of `Y` (pass `NULL` if not needed)
1373: Level: advanced
1375: Notes:
1376: Concatenation is similar to the functionality of a `VECNEST` object; they both represent combination of
1377: different vector spaces. However, concatenated vectors do not store any information about their
1378: sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1379: manipulation of data in the concatenated vector that corresponds to the original components at creation.
1381: This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1382: has to operate on combined vector spaces and cannot utilize `VECNEST` objects due to incompatibility with
1383: bound projections.
1385: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1386: @*/
1387: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1388: {
1389: MPI_Comm comm;
1390: VecType vec_type;
1391: Vec Ytmp, Xtmp;
1392: IS *is_tmp;
1393: PetscInt i, shift = 0, Xnl, Xng, Xbegin;
1395: PetscFunctionBegin;
1399: PetscAssertPointer(Y, 3);
1401: if ((*X)->ops->concatenate) {
1402: /* use the dedicated concatenation function if available */
1403: PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1404: } else {
1405: /* loop over vectors and start creating IS */
1406: comm = PetscObjectComm((PetscObject)*X);
1407: PetscCall(VecGetType(*X, &vec_type));
1408: PetscCall(PetscMalloc1(nx, &is_tmp));
1409: for (i = 0; i < nx; i++) {
1410: PetscCall(VecGetSize(X[i], &Xng));
1411: PetscCall(VecGetLocalSize(X[i], &Xnl));
1412: PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1413: PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1414: shift += Xng;
1415: }
1416: /* create the concatenated vector */
1417: PetscCall(VecCreate(comm, &Ytmp));
1418: PetscCall(VecSetType(Ytmp, vec_type));
1419: PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1420: PetscCall(VecSetUp(Ytmp));
1421: /* copy data from X array to Y and return */
1422: for (i = 0; i < nx; i++) {
1423: PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1424: PetscCall(VecCopy(X[i], Xtmp));
1425: PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1426: }
1427: *Y = Ytmp;
1428: if (x_is) {
1429: *x_is = is_tmp;
1430: } else {
1431: for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1432: PetscCall(PetscFree(is_tmp));
1433: }
1434: }
1435: PetscFunctionReturn(PETSC_SUCCESS);
1436: }
1438: /* A helper function for VecGetSubVector to check if we can implement it with no-copy (i.e. the subvector shares
1439: memory with the original vector), and the block size of the subvector.
1441: Input Parameters:
1442: + X - the original vector
1443: - is - the index set of the subvector
1445: Output Parameters:
1446: + contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
1447: . start - start of contiguous block, as an offset from the start of the ownership range of the original vector
1448: - blocksize - the block size of the subvector
1450: */
1451: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1452: {
1453: PetscInt gstart, gend, lstart;
1454: PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1455: PetscInt n, N, ibs, vbs, bs = -1;
1457: PetscFunctionBegin;
1458: PetscCall(ISGetLocalSize(is, &n));
1459: PetscCall(ISGetSize(is, &N));
1460: PetscCall(ISGetBlockSize(is, &ibs));
1461: PetscCall(VecGetBlockSize(X, &vbs));
1462: PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1463: PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1464: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1465: if (ibs > 1) {
1466: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1467: bs = ibs;
1468: } else {
1469: if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1470: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1471: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1472: }
1474: *contig = red[0];
1475: *start = lstart;
1476: *blocksize = bs;
1477: PetscFunctionReturn(PETSC_SUCCESS);
1478: }
1480: /* A helper function for VecGetSubVector, to be used when we have to build a standalone subvector through VecScatter
1482: Input Parameters:
1483: + X - the original vector
1484: . is - the index set of the subvector
1485: - bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()
1487: Output Parameter:
1488: . Z - the subvector, which will compose the VecScatter context on output
1489: */
1490: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1491: {
1492: PetscInt n, N;
1493: VecScatter vscat;
1494: Vec Y;
1496: PetscFunctionBegin;
1497: PetscCall(ISGetLocalSize(is, &n));
1498: PetscCall(ISGetSize(is, &N));
1499: PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1500: PetscCall(VecSetSizes(Y, n, N));
1501: PetscCall(VecSetBlockSize(Y, bs));
1502: PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1503: PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1504: PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1505: PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1506: PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1507: PetscCall(VecScatterDestroy(&vscat));
1508: *Z = Y;
1509: PetscFunctionReturn(PETSC_SUCCESS);
1510: }
1512: /*@
1513: VecGetSubVector - Gets a vector representing part of another vector
1515: Collective
1517: Input Parameters:
1518: + X - vector from which to extract a subvector
1519: - is - index set representing portion of `X` to extract
1521: Output Parameter:
1522: . Y - subvector corresponding to `is`
1524: Level: advanced
1526: Notes:
1527: The subvector `Y` should be returned with `VecRestoreSubVector()`.
1528: `X` and `is` must be defined on the same communicator
1530: Changes to the subvector will be reflected in the `X` vector on the call to `VecRestoreSubVector()`.
1532: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1533: modifying the subvector. Other non-overlapping subvectors can still be obtained from `X` using this function.
1535: 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`.
1537: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1538: @*/
1539: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1540: {
1541: Vec Z;
1543: PetscFunctionBegin;
1546: PetscCheckSameComm(X, 1, is, 2);
1547: PetscAssertPointer(Y, 3);
1548: if (X->ops->getsubvector) {
1549: PetscUseTypeMethod(X, getsubvector, is, &Z);
1550: } else { /* Default implementation currently does no caching */
1551: PetscBool contig;
1552: PetscInt n, N, start, bs;
1554: PetscCall(ISGetLocalSize(is, &n));
1555: PetscCall(ISGetSize(is, &N));
1556: PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1557: if (contig) { /* We can do a no-copy implementation */
1558: const PetscScalar *x;
1559: PetscInt state = 0;
1560: PetscBool isstd, iscuda, iship;
1562: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1563: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1564: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1565: if (iscuda) {
1566: #if defined(PETSC_HAVE_CUDA)
1567: const PetscScalar *x_d;
1568: PetscMPIInt size;
1569: PetscOffloadMask flg;
1571: PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1572: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1573: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1574: if (x) x += start;
1575: if (x_d) x_d += start;
1576: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1577: if (size == 1) {
1578: PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1579: } else {
1580: PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1581: }
1582: Z->offloadmask = flg;
1583: #endif
1584: } else if (iship) {
1585: #if defined(PETSC_HAVE_HIP)
1586: const PetscScalar *x_d;
1587: PetscMPIInt size;
1588: PetscOffloadMask flg;
1590: PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1591: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1592: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1593: if (x) x += start;
1594: if (x_d) x_d += start;
1595: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1596: if (size == 1) {
1597: PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1598: } else {
1599: PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1600: }
1601: Z->offloadmask = flg;
1602: #endif
1603: } else if (isstd) {
1604: PetscMPIInt size;
1606: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1607: PetscCall(VecGetArrayRead(X, &x));
1608: if (x) x += start;
1609: if (size == 1) {
1610: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1611: } else {
1612: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1613: }
1614: PetscCall(VecRestoreArrayRead(X, &x));
1615: } else { /* default implementation: use place array */
1616: PetscCall(VecGetArrayRead(X, &x));
1617: PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1618: PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1619: PetscCall(VecSetSizes(Z, n, N));
1620: PetscCall(VecSetBlockSize(Z, bs));
1621: PetscCall(VecPlaceArray(Z, PetscSafePointerPlusOffset(x, start)));
1622: PetscCall(VecRestoreArrayRead(X, &x));
1623: }
1625: /* this is relevant only in debug mode */
1626: PetscCall(VecLockGet(X, &state));
1627: if (state) PetscCall(VecLockReadPush(Z));
1628: Z->ops->placearray = NULL;
1629: Z->ops->replacearray = NULL;
1630: } else { /* Have to create a scatter and do a copy */
1631: PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1632: }
1633: }
1634: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1635: if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1636: PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1637: *Y = Z;
1638: PetscFunctionReturn(PETSC_SUCCESS);
1639: }
1641: /*@
1642: VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`
1644: Collective
1646: Input Parameters:
1647: + X - vector from which subvector was obtained
1648: . is - index set representing the subset of `X`
1649: - Y - subvector being restored
1651: Level: advanced
1653: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1654: @*/
1655: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1656: {
1657: PETSC_UNUSED PetscObjectState dummystate = 0;
1658: PetscBool unchanged;
1660: PetscFunctionBegin;
1663: PetscCheckSameComm(X, 1, is, 2);
1664: PetscAssertPointer(Y, 3);
1667: if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1668: else {
1669: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1670: if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1671: VecScatter scatter;
1672: PetscInt state;
1674: PetscCall(VecLockGet(X, &state));
1675: PetscCheck(state == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec X is locked for read-only or read/write access");
1677: PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1678: if (scatter) {
1679: PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1680: PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1681: } else {
1682: PetscBool iscuda, iship;
1683: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1684: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1686: if (iscuda) {
1687: #if defined(PETSC_HAVE_CUDA)
1688: PetscOffloadMask ymask = (*Y)->offloadmask;
1690: /* The offloadmask of X dictates where to move memory
1691: If X GPU data is valid, then move Y data on GPU if needed
1692: Otherwise, move back to the CPU */
1693: switch (X->offloadmask) {
1694: case PETSC_OFFLOAD_BOTH:
1695: if (ymask == PETSC_OFFLOAD_CPU) {
1696: PetscCall(VecCUDAResetArray(*Y));
1697: } else if (ymask == PETSC_OFFLOAD_GPU) {
1698: X->offloadmask = PETSC_OFFLOAD_GPU;
1699: }
1700: break;
1701: case PETSC_OFFLOAD_GPU:
1702: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1703: break;
1704: case PETSC_OFFLOAD_CPU:
1705: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1706: break;
1707: case PETSC_OFFLOAD_UNALLOCATED:
1708: case PETSC_OFFLOAD_KOKKOS:
1709: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1710: }
1711: #endif
1712: } else if (iship) {
1713: #if defined(PETSC_HAVE_HIP)
1714: PetscOffloadMask ymask = (*Y)->offloadmask;
1716: /* The offloadmask of X dictates where to move memory
1717: If X GPU data is valid, then move Y data on GPU if needed
1718: Otherwise, move back to the CPU */
1719: switch (X->offloadmask) {
1720: case PETSC_OFFLOAD_BOTH:
1721: if (ymask == PETSC_OFFLOAD_CPU) {
1722: PetscCall(VecHIPResetArray(*Y));
1723: } else if (ymask == PETSC_OFFLOAD_GPU) {
1724: X->offloadmask = PETSC_OFFLOAD_GPU;
1725: }
1726: break;
1727: case PETSC_OFFLOAD_GPU:
1728: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1729: break;
1730: case PETSC_OFFLOAD_CPU:
1731: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1732: break;
1733: case PETSC_OFFLOAD_UNALLOCATED:
1734: case PETSC_OFFLOAD_KOKKOS:
1735: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1736: }
1737: #endif
1738: } else {
1739: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1740: PetscCall(VecResetArray(*Y));
1741: }
1742: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1743: }
1744: }
1745: }
1746: PetscCall(VecDestroy(Y));
1747: PetscFunctionReturn(PETSC_SUCCESS);
1748: }
1750: /*@
1751: VecCreateLocalVector - Creates a vector object suitable for use with `VecGetLocalVector()` and friends. You must call `VecDestroy()` when the
1752: vector is no longer needed.
1754: Not Collective.
1756: Input Parameter:
1757: . v - The vector for which the local vector is desired.
1759: Output Parameter:
1760: . w - Upon exit this contains the local vector.
1762: Level: beginner
1764: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1765: @*/
1766: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1767: {
1768: PetscMPIInt size;
1770: PetscFunctionBegin;
1772: PetscAssertPointer(w, 2);
1773: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1774: if (size == 1) PetscCall(VecDuplicate(v, w));
1775: else if (v->ops->createlocalvector) PetscUseTypeMethod(v, createlocalvector, w);
1776: else {
1777: VecType type;
1778: PetscInt n;
1780: PetscCall(VecCreate(PETSC_COMM_SELF, w));
1781: PetscCall(VecGetLocalSize(v, &n));
1782: PetscCall(VecSetSizes(*w, n, n));
1783: PetscCall(VecGetBlockSize(v, &n));
1784: PetscCall(VecSetBlockSize(*w, n));
1785: PetscCall(VecGetType(v, &type));
1786: PetscCall(VecSetType(*w, type));
1787: }
1788: PetscFunctionReturn(PETSC_SUCCESS);
1789: }
1791: /*@
1792: VecGetLocalVectorRead - Maps the local portion of a vector into a
1793: vector.
1795: Not Collective.
1797: Input Parameter:
1798: . v - The vector for which the local vector is desired.
1800: Output Parameter:
1801: . w - Upon exit this contains the local vector.
1803: Level: beginner
1805: Notes:
1806: You must call `VecRestoreLocalVectorRead()` when the local
1807: vector is no longer needed.
1809: This function is similar to `VecGetArrayRead()` which maps the local
1810: portion into a raw pointer. `VecGetLocalVectorRead()` is usually
1811: almost as efficient as `VecGetArrayRead()` but in certain circumstances
1812: `VecGetLocalVectorRead()` can be much more efficient than
1813: `VecGetArrayRead()`. This is because the construction of a contiguous
1814: array representing the vector data required by `VecGetArrayRead()` can
1815: be an expensive operation for certain vector types. For example, for
1816: GPU vectors `VecGetArrayRead()` requires that the data between device
1817: and host is synchronized.
1819: Unlike `VecGetLocalVector()`, this routine is not collective and
1820: preserves cached information.
1822: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1823: @*/
1824: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1825: {
1826: PetscFunctionBegin;
1829: VecCheckSameLocalSize(v, 1, w, 2);
1830: if (v->ops->getlocalvectorread) {
1831: PetscUseTypeMethod(v, getlocalvectorread, w);
1832: } else {
1833: PetscScalar *a;
1835: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1836: PetscCall(VecPlaceArray(w, a));
1837: }
1838: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1839: PetscCall(VecLockReadPush(v));
1840: PetscCall(VecLockReadPush(w));
1841: PetscFunctionReturn(PETSC_SUCCESS);
1842: }
1844: /*@
1845: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1846: previously mapped into a vector using `VecGetLocalVectorRead()`.
1848: Not Collective.
1850: Input Parameters:
1851: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
1852: - w - The vector into which the local portion of `v` was mapped.
1854: Level: beginner
1856: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1857: @*/
1858: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1859: {
1860: PetscFunctionBegin;
1863: if (v->ops->restorelocalvectorread) {
1864: PetscUseTypeMethod(v, restorelocalvectorread, w);
1865: } else {
1866: const PetscScalar *a;
1868: PetscCall(VecGetArrayRead(w, &a));
1869: PetscCall(VecRestoreArrayRead(v, &a));
1870: PetscCall(VecResetArray(w));
1871: }
1872: PetscCall(VecLockReadPop(v));
1873: PetscCall(VecLockReadPop(w));
1874: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1875: PetscFunctionReturn(PETSC_SUCCESS);
1876: }
1878: /*@
1879: VecGetLocalVector - Maps the local portion of a vector into a
1880: vector.
1882: Collective
1884: Input Parameter:
1885: . v - The vector for which the local vector is desired.
1887: Output Parameter:
1888: . w - Upon exit this contains the local vector.
1890: Level: beginner
1892: Notes:
1893: You must call `VecRestoreLocalVector()` when the local
1894: vector is no longer needed.
1896: This function is similar to `VecGetArray()` which maps the local
1897: portion into a raw pointer. `VecGetLocalVector()` is usually about as
1898: efficient as `VecGetArray()` but in certain circumstances
1899: `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1900: This is because the construction of a contiguous array representing
1901: the vector data required by `VecGetArray()` can be an expensive
1902: operation for certain vector types. For example, for GPU vectors
1903: `VecGetArray()` requires that the data between device and host is
1904: synchronized.
1906: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1907: @*/
1908: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1909: {
1910: PetscFunctionBegin;
1913: VecCheckSameLocalSize(v, 1, w, 2);
1914: if (v->ops->getlocalvector) {
1915: PetscUseTypeMethod(v, getlocalvector, w);
1916: } else {
1917: PetscScalar *a;
1919: PetscCall(VecGetArray(v, &a));
1920: PetscCall(VecPlaceArray(w, a));
1921: }
1922: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1923: PetscFunctionReturn(PETSC_SUCCESS);
1924: }
1926: /*@
1927: VecRestoreLocalVector - Unmaps the local portion of a vector
1928: previously mapped into a vector using `VecGetLocalVector()`.
1930: Logically Collective.
1932: Input Parameters:
1933: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
1934: - w - The vector into which the local portion of `v` was mapped.
1936: Level: beginner
1938: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1939: @*/
1940: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1941: {
1942: PetscFunctionBegin;
1945: if (v->ops->restorelocalvector) {
1946: PetscUseTypeMethod(v, restorelocalvector, w);
1947: } else {
1948: PetscScalar *a;
1949: PetscCall(VecGetArray(w, &a));
1950: PetscCall(VecRestoreArray(v, &a));
1951: PetscCall(VecResetArray(w));
1952: }
1953: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1954: PetscCall(PetscObjectStateIncrease((PetscObject)v));
1955: PetscFunctionReturn(PETSC_SUCCESS);
1956: }
1958: /*@C
1959: VecGetArray - Returns a pointer to a contiguous array that contains this
1960: MPI processes's portion of the vector data
1962: Logically Collective
1964: Input Parameter:
1965: . x - the vector
1967: Output Parameter:
1968: . a - location to put pointer to the array
1970: Level: beginner
1972: Notes:
1973: For the standard PETSc vectors, `VecGetArray()` returns a pointer to the local data array and
1974: does not use any copies. If the underlying vector data is not stored in a contiguous array
1975: this routine will copy the data to a contiguous array and return a pointer to that. You MUST
1976: call `VecRestoreArray()` when you no longer need access to the array.
1978: Fortran Notes:
1979: `VecGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayF90()`
1981: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
1982: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
1983: @*/
1984: PetscErrorCode VecGetArray(Vec x, PetscScalar **a)
1985: {
1986: PetscFunctionBegin;
1988: PetscCall(VecSetErrorIfLocked(x, 1));
1989: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
1990: PetscUseTypeMethod(x, getarray, a);
1991: } else if (x->petscnative) { /* VECSTANDARD */
1992: *a = *((PetscScalar **)x->data);
1993: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
1994: PetscFunctionReturn(PETSC_SUCCESS);
1995: }
1997: /*@C
1998: VecRestoreArray - Restores a vector after `VecGetArray()` has been called and the array is no longer needed
2000: Logically Collective
2002: Input Parameters:
2003: + x - the vector
2004: - a - location of pointer to array obtained from `VecGetArray()`
2006: Level: beginner
2008: Fortran Notes:
2009: `VecRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayF90()`
2011: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2012: `VecGetArrayPair()`, `VecRestoreArrayPair()`
2013: @*/
2014: PetscErrorCode VecRestoreArray(Vec x, PetscScalar **a)
2015: {
2016: PetscFunctionBegin;
2018: if (a) PetscAssertPointer(a, 2);
2019: if (x->ops->restorearray) {
2020: PetscUseTypeMethod(x, restorearray, a);
2021: } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
2022: if (a) *a = NULL;
2023: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2024: PetscFunctionReturn(PETSC_SUCCESS);
2025: }
2026: /*@C
2027: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
2029: Not Collective
2031: Input Parameter:
2032: . x - the vector
2034: Output Parameter:
2035: . a - the array
2037: Level: beginner
2039: Notes:
2040: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2042: Unlike `VecGetArray()`, preserves cached information like vector norms.
2044: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
2045: implementations may require a copy, but such implementations should cache the contiguous representation so that
2046: only one copy is performed when this routine is called multiple times in sequence.
2048: Fortran Notes:
2049: `VecGetArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayReadF90()`
2051: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2052: @*/
2053: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar **a)
2054: {
2055: PetscFunctionBegin;
2057: PetscAssertPointer(a, 2);
2058: if (x->ops->getarrayread) {
2059: PetscUseTypeMethod(x, getarrayread, a);
2060: } else if (x->ops->getarray) {
2061: PetscObjectState state;
2063: /* VECNEST, VECCUDA, VECKOKKOS etc */
2064: // x->ops->getarray may bump the object state, but since we know this is a read-only get
2065: // we can just undo that
2066: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2067: PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
2068: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2069: } else if (x->petscnative) {
2070: /* VECSTANDARD */
2071: *a = *((PetscScalar **)x->data);
2072: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
2073: PetscFunctionReturn(PETSC_SUCCESS);
2074: }
2076: /*@C
2077: VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`
2079: Not Collective
2081: Input Parameters:
2082: + x - the vector
2083: - a - the array
2085: Level: beginner
2087: Fortran Notes:
2088: `VecRestoreArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayReadF90()`
2090: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2091: @*/
2092: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar **a)
2093: {
2094: PetscFunctionBegin;
2096: if (a) PetscAssertPointer(a, 2);
2097: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
2098: /* nothing */
2099: } else if (x->ops->restorearrayread) { /* VECNEST */
2100: PetscUseTypeMethod(x, restorearrayread, a);
2101: } else { /* No one? */
2102: PetscObjectState state;
2104: // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2105: // we can just undo that
2106: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2107: PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2108: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2109: }
2110: if (a) *a = NULL;
2111: PetscFunctionReturn(PETSC_SUCCESS);
2112: }
2114: /*@C
2115: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
2116: MPI processes's portion of the vector data.
2118: Logically Collective
2120: Input Parameter:
2121: . x - the vector
2123: Output Parameter:
2124: . a - location to put pointer to the array
2126: Level: intermediate
2128: Note:
2129: The values in this array are NOT valid, the caller of this routine is responsible for putting
2130: values into the array; any values it does not set will be invalid.
2132: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2134: For vectors associated with GPUs, the host and device vectors are not synchronized before
2135: giving access. If you need correct values in the array use `VecGetArray()`
2137: Fortran Notes:
2138: `VecGetArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayWriteF90()`
2140: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteF90()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
2141: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
2142: @*/
2143: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar **a)
2144: {
2145: PetscFunctionBegin;
2147: PetscAssertPointer(a, 2);
2148: PetscCall(VecSetErrorIfLocked(x, 1));
2149: if (x->ops->getarraywrite) {
2150: PetscUseTypeMethod(x, getarraywrite, a);
2151: } else {
2152: PetscCall(VecGetArray(x, a));
2153: }
2154: PetscFunctionReturn(PETSC_SUCCESS);
2155: }
2157: /*@C
2158: VecRestoreArrayWrite - Restores a vector after `VecGetArrayWrite()` has been called.
2160: Logically Collective
2162: Input Parameters:
2163: + x - the vector
2164: - a - location of pointer to array obtained from `VecGetArray()`
2166: Level: beginner
2168: Fortran Notes:
2169: `VecRestoreArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayWriteF90()`
2171: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteF90()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2172: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2173: @*/
2174: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar **a)
2175: {
2176: PetscFunctionBegin;
2178: if (a) PetscAssertPointer(a, 2);
2179: if (x->ops->restorearraywrite) {
2180: PetscUseTypeMethod(x, restorearraywrite, a);
2181: } else if (x->ops->restorearray) {
2182: PetscUseTypeMethod(x, restorearray, a);
2183: }
2184: if (a) *a = NULL;
2185: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2186: PetscFunctionReturn(PETSC_SUCCESS);
2187: }
2189: /*@C
2190: VecGetArrays - Returns a pointer to the arrays in a set of vectors
2191: that were created by a call to `VecDuplicateVecs()`.
2193: Logically Collective; No Fortran Support
2195: Input Parameters:
2196: + x - the vectors
2197: - n - the number of vectors
2199: Output Parameter:
2200: . a - location to put pointer to the array
2202: Level: intermediate
2204: Note:
2205: You MUST call `VecRestoreArrays()` when you no longer need access to the arrays.
2207: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2208: @*/
2209: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2210: {
2211: PetscInt i;
2212: PetscScalar **q;
2214: PetscFunctionBegin;
2215: PetscAssertPointer(x, 1);
2217: PetscAssertPointer(a, 3);
2218: PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2219: PetscCall(PetscMalloc1(n, &q));
2220: for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2221: *a = q;
2222: PetscFunctionReturn(PETSC_SUCCESS);
2223: }
2225: /*@C
2226: VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2227: has been called.
2229: Logically Collective; No Fortran Support
2231: Input Parameters:
2232: + x - the vector
2233: . n - the number of vectors
2234: - a - location of pointer to arrays obtained from `VecGetArrays()`
2236: Notes:
2237: For regular PETSc vectors this routine does not involve any copies. For
2238: any special vectors that do not store local vector data in a contiguous
2239: array, this routine will copy the data back into the underlying
2240: vector data structure from the arrays obtained with `VecGetArrays()`.
2242: Level: intermediate
2244: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2245: @*/
2246: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2247: {
2248: PetscInt i;
2249: PetscScalar **q = *a;
2251: PetscFunctionBegin;
2252: PetscAssertPointer(x, 1);
2254: PetscAssertPointer(a, 3);
2256: for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2257: PetscCall(PetscFree(q));
2258: PetscFunctionReturn(PETSC_SUCCESS);
2259: }
2261: /*@C
2262: VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g.,
2263: `VECCUDA`), the returned pointer will be a device pointer to the device memory that contains
2264: this MPI processes's portion of the vector data.
2266: Logically Collective; No Fortran Support
2268: Input Parameter:
2269: . x - the vector
2271: Output Parameters:
2272: + a - location to put pointer to the array
2273: - mtype - memory type of the array
2275: Level: beginner
2277: Note:
2278: Device data is guaranteed to have the latest value. Otherwise, when this is a host vector
2279: (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host
2280: pointer.
2282: For `VECKOKKOS`, if Kokkos is configured without device (e.g., use serial or openmp), per
2283: this function, the vector works like `VECSEQ`/`VECMPI`; otherwise, it works like `VECCUDA` or
2284: `VECHIP` etc.
2286: Use `VecRestoreArrayAndMemType()` when the array access is no longer needed.
2288: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`,
2289: `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2290: @*/
2291: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2292: {
2293: PetscFunctionBegin;
2296: PetscAssertPointer(a, 2);
2297: if (mtype) PetscAssertPointer(mtype, 3);
2298: PetscCall(VecSetErrorIfLocked(x, 1));
2299: if (x->ops->getarrayandmemtype) {
2300: /* VECCUDA, VECKOKKOS etc */
2301: PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2302: } else {
2303: /* VECSTANDARD, VECNEST, VECVIENNACL */
2304: PetscCall(VecGetArray(x, a));
2305: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2306: }
2307: PetscFunctionReturn(PETSC_SUCCESS);
2308: }
2310: /*@C
2311: VecRestoreArrayAndMemType - Restores a vector after `VecGetArrayAndMemType()` has been called.
2313: Logically Collective; No Fortran Support
2315: Input Parameters:
2316: + x - the vector
2317: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`
2319: Level: beginner
2321: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`,
2322: `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2323: @*/
2324: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar **a)
2325: {
2326: PetscFunctionBegin;
2329: if (a) PetscAssertPointer(a, 2);
2330: if (x->ops->restorearrayandmemtype) {
2331: /* VECCUDA, VECKOKKOS etc */
2332: PetscUseTypeMethod(x, restorearrayandmemtype, a);
2333: } else {
2334: /* VECNEST, VECVIENNACL */
2335: PetscCall(VecRestoreArray(x, a));
2336: } /* VECSTANDARD does nothing */
2337: if (a) *a = NULL;
2338: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2339: PetscFunctionReturn(PETSC_SUCCESS);
2340: }
2342: /*@C
2343: VecGetArrayReadAndMemType - Like `VecGetArrayRead()`, but if the input vector is a device vector, it will return a read-only device pointer.
2344: The returned pointer is guaranteed to point to up-to-date data. For host vectors, it functions as `VecGetArrayRead()`.
2346: Not Collective; No Fortran Support
2348: Input Parameter:
2349: . x - the vector
2351: Output Parameters:
2352: + a - the array
2353: - mtype - memory type of the array
2355: Level: beginner
2357: Notes:
2358: The array must be returned using a matching call to `VecRestoreArrayReadAndMemType()`.
2360: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2361: @*/
2362: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar **a, PetscMemType *mtype)
2363: {
2364: PetscFunctionBegin;
2367: PetscAssertPointer(a, 2);
2368: if (mtype) PetscAssertPointer(mtype, 3);
2369: if (x->ops->getarrayreadandmemtype) {
2370: /* VECCUDA/VECHIP though they are also petscnative */
2371: PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2372: } else if (x->ops->getarrayandmemtype) {
2373: /* VECKOKKOS */
2374: PetscObjectState state;
2376: // see VecGetArrayRead() for why
2377: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2378: PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2379: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2380: } else {
2381: PetscCall(VecGetArrayRead(x, a));
2382: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2383: }
2384: PetscFunctionReturn(PETSC_SUCCESS);
2385: }
2387: /*@C
2388: VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`
2390: Not Collective; No Fortran Support
2392: Input Parameters:
2393: + x - the vector
2394: - a - the array
2396: Level: beginner
2398: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2399: @*/
2400: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar **a)
2401: {
2402: PetscFunctionBegin;
2405: if (a) PetscAssertPointer(a, 2);
2406: if (x->ops->restorearrayreadandmemtype) {
2407: /* VECCUDA/VECHIP */
2408: PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2409: } else if (!x->petscnative) {
2410: /* VECNEST */
2411: PetscCall(VecRestoreArrayRead(x, a));
2412: }
2413: if (a) *a = NULL;
2414: PetscFunctionReturn(PETSC_SUCCESS);
2415: }
2417: /*@C
2418: VecGetArrayWriteAndMemType - Like `VecGetArrayWrite()`, but if this is a device vector it will always return
2419: a device pointer to the device memory that contains this processor's portion of the vector data.
2421: Logically Collective; No Fortran Support
2423: Input Parameter:
2424: . x - the vector
2426: Output Parameters:
2427: + a - the array
2428: - mtype - memory type of the array
2430: Level: beginner
2432: Note:
2433: The array must be returned using a matching call to `VecRestoreArrayWriteAndMemType()`, where it will label the device memory as most recent.
2435: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2436: @*/
2437: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2438: {
2439: PetscFunctionBegin;
2442: PetscCall(VecSetErrorIfLocked(x, 1));
2443: PetscAssertPointer(a, 2);
2444: if (mtype) PetscAssertPointer(mtype, 3);
2445: if (x->ops->getarraywriteandmemtype) {
2446: /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2447: PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2448: } else if (x->ops->getarrayandmemtype) {
2449: PetscCall(VecGetArrayAndMemType(x, a, mtype));
2450: } else {
2451: /* VECNEST, VECVIENNACL */
2452: PetscCall(VecGetArrayWrite(x, a));
2453: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2454: }
2455: PetscFunctionReturn(PETSC_SUCCESS);
2456: }
2458: /*@C
2459: VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`
2461: Logically Collective; No Fortran Support
2463: Input Parameters:
2464: + x - the vector
2465: - a - the array
2467: Level: beginner
2469: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2470: @*/
2471: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar **a)
2472: {
2473: PetscFunctionBegin;
2476: PetscCall(VecSetErrorIfLocked(x, 1));
2477: if (a) PetscAssertPointer(a, 2);
2478: if (x->ops->restorearraywriteandmemtype) {
2479: /* VECCUDA/VECHIP */
2480: PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2481: PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2482: } else if (x->ops->restorearrayandmemtype) {
2483: PetscCall(VecRestoreArrayAndMemType(x, a));
2484: } else {
2485: PetscCall(VecRestoreArray(x, a));
2486: }
2487: if (a) *a = NULL;
2488: PetscFunctionReturn(PETSC_SUCCESS);
2489: }
2491: /*@
2492: VecPlaceArray - Allows one to replace the array in a vector with an
2493: array provided by the user. This is useful to avoid copying an array
2494: into a vector.
2496: Logically Collective; No Fortran Support
2498: Input Parameters:
2499: + vec - the vector
2500: - array - the array
2502: Level: developer
2504: Notes:
2505: Use `VecReplaceArray()` instead to permanently replace the array
2507: You can return to the original array with a call to `VecResetArray()`. `vec` does not take
2508: ownership of `array` in any way.
2510: The user must free `array` themselves but be careful not to
2511: do so before the vector has either been destroyed, had its original array restored with
2512: `VecResetArray()` or permanently replaced with `VecReplaceArray()`.
2514: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2515: @*/
2516: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2517: {
2518: PetscFunctionBegin;
2521: if (array) PetscAssertPointer(array, 2);
2522: PetscUseTypeMethod(vec, placearray, array);
2523: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2524: PetscFunctionReturn(PETSC_SUCCESS);
2525: }
2527: /*@C
2528: VecReplaceArray - Allows one to replace the array in a vector with an
2529: array provided by the user. This is useful to avoid copying an array
2530: into a vector.
2532: Logically Collective; No Fortran Support
2534: Input Parameters:
2535: + vec - the vector
2536: - array - the array
2538: Level: developer
2540: Notes:
2541: This permanently replaces the array and frees the memory associated
2542: with the old array. Use `VecPlaceArray()` to temporarily replace the array.
2544: The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
2545: freed by the user. It will be freed when the vector is destroyed.
2547: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2548: @*/
2549: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2550: {
2551: PetscFunctionBegin;
2554: PetscUseTypeMethod(vec, replacearray, array);
2555: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2556: PetscFunctionReturn(PETSC_SUCCESS);
2557: }
2559: /*MC
2560: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
2561: and makes them accessible via a Fortran pointer.
2563: Synopsis:
2564: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
2566: Collective
2568: Input Parameters:
2569: + x - a vector to mimic
2570: - n - the number of vectors to obtain
2572: Output Parameters:
2573: + y - Fortran pointer to the array of vectors
2574: - ierr - error code
2576: Example of Usage:
2577: .vb
2578: #include <petsc/finclude/petscvec.h>
2579: use petscvec
2581: Vec x
2582: Vec, pointer :: y(:)
2583: ....
2584: call VecDuplicateVecsF90(x,2,y,ierr)
2585: call VecSet(y(2),alpha,ierr)
2586: call VecSet(y(2),alpha,ierr)
2587: ....
2588: call VecDestroyVecsF90(2,y,ierr)
2589: .ve
2591: Level: beginner
2593: Note:
2594: Use `VecDestroyVecsF90()` to free the space.
2596: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecsF90()`, `VecDuplicateVecs()`
2597: M*/
2599: /*MC
2600: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2601: `VecGetArrayF90()`.
2603: Synopsis:
2604: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2606: Logically Collective
2608: Input Parameters:
2609: + x - vector
2610: - xx_v - the Fortran pointer to the array
2612: Output Parameter:
2613: . ierr - error code
2615: Example of Usage:
2616: .vb
2617: #include <petsc/finclude/petscvec.h>
2618: use petscvec
2620: PetscScalar, pointer :: xx_v(:)
2621: ....
2622: call VecGetArrayF90(x,xx_v,ierr)
2623: xx_v(3) = a
2624: call VecRestoreArrayF90(x,xx_v,ierr)
2625: .ve
2627: Level: beginner
2629: .seealso: [](ch_vectors), `Vec`, `VecGetArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrayReadF90()`
2630: M*/
2632: /*MC
2633: VecDestroyVecsF90 - Frees a block of vectors obtained with `VecDuplicateVecsF90()`.
2635: Synopsis:
2636: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
2638: Collective
2640: Input Parameters:
2641: + n - the number of vectors previously obtained
2642: - x - pointer to array of vector pointers
2644: Output Parameter:
2645: . ierr - error code
2647: Level: beginner
2649: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecs()`, `VecDuplicateVecsF90()`
2650: M*/
2652: /*MC
2653: VecGetArrayF90 - Accesses a vector array from Fortran. For default PETSc
2654: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2655: this routine is implementation dependent. You MUST call `VecRestoreArrayF90()`
2656: when you no longer need access to the array.
2658: Synopsis:
2659: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2661: Logically Collective
2663: Input Parameter:
2664: . x - vector
2666: Output Parameters:
2667: + xx_v - the Fortran pointer to the array
2668: - ierr - error code
2670: Example of Usage:
2671: .vb
2672: #include <petsc/finclude/petscvec.h>
2673: use petscvec
2675: PetscScalar, pointer :: xx_v(:)
2676: ....
2677: call VecGetArrayF90(x,xx_v,ierr)
2678: xx_v(3) = a
2679: call VecRestoreArrayF90(x,xx_v,ierr)
2680: .ve
2682: Level: beginner
2684: Note:
2685: If you ONLY intend to read entries from the array and not change any entries you should use `VecGetArrayReadF90()`.
2687: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayReadF90()`
2688: M*/
2690: /*MC
2691: VecGetArrayReadF90 - Accesses a read only array from Fortran. For default PETSc
2692: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2693: this routine is implementation dependent. You MUST call `VecRestoreArrayReadF90()`
2694: when you no longer need access to the array.
2696: Synopsis:
2697: VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2699: Logically Collective
2701: Input Parameter:
2702: . x - vector
2704: Output Parameters:
2705: + xx_v - the Fortran pointer to the array
2706: - ierr - error code
2708: Example of Usage:
2709: .vb
2710: #include <petsc/finclude/petscvec.h>
2711: use petscvec
2713: PetscScalar, pointer :: xx_v(:)
2714: ....
2715: call VecGetArrayReadF90(x,xx_v,ierr)
2716: a = xx_v(3)
2717: call VecRestoreArrayReadF90(x,xx_v,ierr)
2718: .ve
2720: Level: beginner
2722: Note:
2723: If you intend to write entries into the array you must use `VecGetArrayF90()`.
2725: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecGetArrayF90()`
2726: M*/
2728: /*MC
2729: VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2730: `VecGetArrayReadF90()`.
2732: Synopsis:
2733: VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2735: Logically Collective
2737: Input Parameters:
2738: + x - vector
2739: - xx_v - the Fortran pointer to the array
2741: Output Parameter:
2742: . ierr - error code
2744: Example of Usage:
2745: .vb
2746: #include <petsc/finclude/petscvec.h>
2747: use petscvec
2749: PetscScalar, pointer :: xx_v(:)
2750: ....
2751: call VecGetArrayReadF90(x,xx_v,ierr)
2752: a = xx_v(3)
2753: call VecRestoreArrayReadF90(x,xx_v,ierr)
2754: .ve
2756: Level: beginner
2758: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecRestoreArrayF90()`
2759: M*/
2761: /*@C
2762: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2763: processor's portion of the vector data. You MUST call `VecRestoreArray2d()`
2764: when you no longer need access to the array.
2766: Logically Collective
2768: Input Parameters:
2769: + x - the vector
2770: . m - first dimension of two dimensional array
2771: . n - second dimension of two dimensional array
2772: . mstart - first index you will use in first coordinate direction (often 0)
2773: - nstart - first index in the second coordinate direction (often 0)
2775: Output Parameter:
2776: . a - location to put pointer to the array
2778: Level: developer
2780: Notes:
2781: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2782: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2783: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2784: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2786: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2788: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2789: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2790: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2791: @*/
2792: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2793: {
2794: PetscInt i, N;
2795: PetscScalar *aa;
2797: PetscFunctionBegin;
2799: PetscAssertPointer(a, 6);
2801: PetscCall(VecGetLocalSize(x, &N));
2802: 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);
2803: PetscCall(VecGetArray(x, &aa));
2805: PetscCall(PetscMalloc1(m, a));
2806: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2807: *a -= mstart;
2808: PetscFunctionReturn(PETSC_SUCCESS);
2809: }
2811: /*@C
2812: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2813: processor's portion of the vector data. You MUST call `VecRestoreArray2dWrite()`
2814: when you no longer need access to the array.
2816: Logically Collective
2818: Input Parameters:
2819: + x - the vector
2820: . m - first dimension of two dimensional array
2821: . n - second dimension of two dimensional array
2822: . mstart - first index you will use in first coordinate direction (often 0)
2823: - nstart - first index in the second coordinate direction (often 0)
2825: Output Parameter:
2826: . a - location to put pointer to the array
2828: Level: developer
2830: Notes:
2831: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2832: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2833: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2834: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2836: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2838: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2839: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2840: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2841: @*/
2842: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2843: {
2844: PetscInt i, N;
2845: PetscScalar *aa;
2847: PetscFunctionBegin;
2849: PetscAssertPointer(a, 6);
2851: PetscCall(VecGetLocalSize(x, &N));
2852: 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);
2853: PetscCall(VecGetArrayWrite(x, &aa));
2855: PetscCall(PetscMalloc1(m, a));
2856: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2857: *a -= mstart;
2858: PetscFunctionReturn(PETSC_SUCCESS);
2859: }
2861: /*@C
2862: VecRestoreArray2d - Restores a vector after `VecGetArray2d()` has been called.
2864: Logically Collective
2866: Input Parameters:
2867: + x - the vector
2868: . m - first dimension of two dimensional array
2869: . n - second dimension of the two dimensional array
2870: . mstart - first index you will use in first coordinate direction (often 0)
2871: . nstart - first index in the second coordinate direction (often 0)
2872: - a - location of pointer to array obtained from `VecGetArray2d()`
2874: Level: developer
2876: Notes:
2877: For regular PETSc vectors this routine does not involve any copies. For
2878: any special vectors that do not store local vector data in a contiguous
2879: array, this routine will copy the data back into the underlying
2880: vector data structure from the array obtained with `VecGetArray()`.
2882: This routine actually zeros out the `a` pointer.
2884: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2885: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2886: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2887: @*/
2888: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2889: {
2890: void *dummy;
2892: PetscFunctionBegin;
2894: PetscAssertPointer(a, 6);
2896: dummy = (void *)(*a + mstart);
2897: PetscCall(PetscFree(dummy));
2898: PetscCall(VecRestoreArray(x, NULL));
2899: *a = NULL;
2900: PetscFunctionReturn(PETSC_SUCCESS);
2901: }
2903: /*@C
2904: VecRestoreArray2dWrite - Restores a vector after `VecGetArray2dWrite()` has been called.
2906: Logically Collective
2908: Input Parameters:
2909: + x - the vector
2910: . m - first dimension of two dimensional array
2911: . n - second dimension of the two dimensional array
2912: . mstart - first index you will use in first coordinate direction (often 0)
2913: . nstart - first index in the second coordinate direction (often 0)
2914: - a - location of pointer to array obtained from `VecGetArray2d()`
2916: Level: developer
2918: Notes:
2919: For regular PETSc vectors this routine does not involve any copies. For
2920: any special vectors that do not store local vector data in a contiguous
2921: array, this routine will copy the data back into the underlying
2922: vector data structure from the array obtained with `VecGetArray()`.
2924: This routine actually zeros out the `a` pointer.
2926: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2927: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2928: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2929: @*/
2930: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2931: {
2932: void *dummy;
2934: PetscFunctionBegin;
2936: PetscAssertPointer(a, 6);
2938: dummy = (void *)(*a + mstart);
2939: PetscCall(PetscFree(dummy));
2940: PetscCall(VecRestoreArrayWrite(x, NULL));
2941: PetscFunctionReturn(PETSC_SUCCESS);
2942: }
2944: /*@C
2945: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2946: processor's portion of the vector data. You MUST call `VecRestoreArray1d()`
2947: when you no longer need access to the array.
2949: Logically Collective
2951: Input Parameters:
2952: + x - the vector
2953: . m - first dimension of two dimensional array
2954: - mstart - first index you will use in first coordinate direction (often 0)
2956: Output Parameter:
2957: . a - location to put pointer to the array
2959: Level: developer
2961: Notes:
2962: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2963: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2964: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2966: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2968: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2969: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2970: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2971: @*/
2972: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2973: {
2974: PetscInt N;
2976: PetscFunctionBegin;
2978: PetscAssertPointer(a, 4);
2980: PetscCall(VecGetLocalSize(x, &N));
2981: 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);
2982: PetscCall(VecGetArray(x, a));
2983: *a -= mstart;
2984: PetscFunctionReturn(PETSC_SUCCESS);
2985: }
2987: /*@C
2988: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
2989: processor's portion of the vector data. You MUST call `VecRestoreArray1dWrite()`
2990: when you no longer need access to the array.
2992: Logically Collective
2994: Input Parameters:
2995: + x - the vector
2996: . m - first dimension of two dimensional array
2997: - mstart - first index you will use in first coordinate direction (often 0)
2999: Output Parameter:
3000: . a - location to put pointer to the array
3002: Level: developer
3004: Notes:
3005: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3006: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3007: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3009: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3011: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3012: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3013: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3014: @*/
3015: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3016: {
3017: PetscInt N;
3019: PetscFunctionBegin;
3021: PetscAssertPointer(a, 4);
3023: PetscCall(VecGetLocalSize(x, &N));
3024: 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);
3025: PetscCall(VecGetArrayWrite(x, a));
3026: *a -= mstart;
3027: PetscFunctionReturn(PETSC_SUCCESS);
3028: }
3030: /*@C
3031: VecRestoreArray1d - Restores a vector after `VecGetArray1d()` has been called.
3033: Logically Collective
3035: Input Parameters:
3036: + x - the vector
3037: . m - first dimension of two dimensional array
3038: . mstart - first index you will use in first coordinate direction (often 0)
3039: - a - location of pointer to array obtained from `VecGetArray1d()`
3041: Level: developer
3043: Notes:
3044: For regular PETSc vectors this routine does not involve any copies. For
3045: any special vectors that do not store local vector data in a contiguous
3046: array, this routine will copy the data back into the underlying
3047: vector data structure from the array obtained with `VecGetArray1d()`.
3049: This routine actually zeros out the `a` pointer.
3051: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3052: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3053: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3054: @*/
3055: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3056: {
3057: PetscFunctionBegin;
3060: PetscCall(VecRestoreArray(x, NULL));
3061: *a = NULL;
3062: PetscFunctionReturn(PETSC_SUCCESS);
3063: }
3065: /*@C
3066: VecRestoreArray1dWrite - Restores a vector after `VecGetArray1dWrite()` has been called.
3068: Logically Collective
3070: Input Parameters:
3071: + x - the vector
3072: . m - first dimension of two dimensional array
3073: . mstart - first index you will use in first coordinate direction (often 0)
3074: - a - location of pointer to array obtained from `VecGetArray1d()`
3076: Level: developer
3078: Notes:
3079: For regular PETSc vectors this routine does not involve any copies. For
3080: any special vectors that do not store local vector data in a contiguous
3081: array, this routine will copy the data back into the underlying
3082: vector data structure from the array obtained with `VecGetArray1d()`.
3084: This routine actually zeros out the `a` pointer.
3086: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3087: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3088: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3089: @*/
3090: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3091: {
3092: PetscFunctionBegin;
3095: PetscCall(VecRestoreArrayWrite(x, NULL));
3096: *a = NULL;
3097: PetscFunctionReturn(PETSC_SUCCESS);
3098: }
3100: /*@C
3101: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
3102: processor's portion of the vector data. You MUST call `VecRestoreArray3d()`
3103: when you no longer need access to the array.
3105: Logically Collective
3107: Input Parameters:
3108: + x - the vector
3109: . m - first dimension of three dimensional array
3110: . n - second dimension of three dimensional array
3111: . p - third dimension of three dimensional array
3112: . mstart - first index you will use in first coordinate direction (often 0)
3113: . nstart - first index in the second coordinate direction (often 0)
3114: - pstart - first index in the third coordinate direction (often 0)
3116: Output Parameter:
3117: . a - location to put pointer to the array
3119: Level: developer
3121: Notes:
3122: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3123: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3124: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3125: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3127: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3129: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3130: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
3131: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3132: @*/
3133: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3134: {
3135: PetscInt i, N, j;
3136: PetscScalar *aa, **b;
3138: PetscFunctionBegin;
3140: PetscAssertPointer(a, 8);
3142: PetscCall(VecGetLocalSize(x, &N));
3143: 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);
3144: PetscCall(VecGetArray(x, &aa));
3146: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3147: b = (PetscScalar **)((*a) + m);
3148: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3149: for (i = 0; i < m; i++)
3150: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset(aa, i * n * p + j * p - pstart);
3151: *a -= mstart;
3152: PetscFunctionReturn(PETSC_SUCCESS);
3153: }
3155: /*@C
3156: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
3157: processor's portion of the vector data. You MUST call `VecRestoreArray3dWrite()`
3158: when you no longer need access to the array.
3160: Logically Collective
3162: Input Parameters:
3163: + x - the vector
3164: . m - first dimension of three dimensional array
3165: . n - second dimension of three dimensional array
3166: . p - third dimension of three dimensional array
3167: . mstart - first index you will use in first coordinate direction (often 0)
3168: . nstart - first index in the second coordinate direction (often 0)
3169: - pstart - first index in the third coordinate direction (often 0)
3171: Output Parameter:
3172: . a - location to put pointer to the array
3174: Level: developer
3176: Notes:
3177: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3178: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3179: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3180: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3182: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3184: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3185: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3186: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3187: @*/
3188: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3189: {
3190: PetscInt i, N, j;
3191: PetscScalar *aa, **b;
3193: PetscFunctionBegin;
3195: PetscAssertPointer(a, 8);
3197: PetscCall(VecGetLocalSize(x, &N));
3198: 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);
3199: PetscCall(VecGetArrayWrite(x, &aa));
3201: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3202: b = (PetscScalar **)((*a) + m);
3203: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3204: for (i = 0; i < m; i++)
3205: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3207: *a -= mstart;
3208: PetscFunctionReturn(PETSC_SUCCESS);
3209: }
3211: /*@C
3212: VecRestoreArray3d - Restores a vector after `VecGetArray3d()` has been called.
3214: Logically Collective
3216: Input Parameters:
3217: + x - the vector
3218: . m - first dimension of three dimensional array
3219: . n - second dimension of the three dimensional array
3220: . p - third dimension of the three dimensional array
3221: . mstart - first index you will use in first coordinate direction (often 0)
3222: . nstart - first index in the second coordinate direction (often 0)
3223: . pstart - first index in the third coordinate direction (often 0)
3224: - a - location of pointer to array obtained from VecGetArray3d()
3226: Level: developer
3228: Notes:
3229: For regular PETSc vectors this routine does not involve any copies. For
3230: any special vectors that do not store local vector data in a contiguous
3231: array, this routine will copy the data back into the underlying
3232: vector data structure from the array obtained with `VecGetArray()`.
3234: This routine actually zeros out the `a` pointer.
3236: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3237: `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3238: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3239: @*/
3240: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3241: {
3242: void *dummy;
3244: PetscFunctionBegin;
3246: PetscAssertPointer(a, 8);
3248: dummy = (void *)(*a + mstart);
3249: PetscCall(PetscFree(dummy));
3250: PetscCall(VecRestoreArray(x, NULL));
3251: *a = NULL;
3252: PetscFunctionReturn(PETSC_SUCCESS);
3253: }
3255: /*@C
3256: VecRestoreArray3dWrite - Restores a vector after `VecGetArray3dWrite()` has been called.
3258: Logically Collective
3260: Input Parameters:
3261: + x - the vector
3262: . m - first dimension of three dimensional array
3263: . n - second dimension of the three dimensional array
3264: . p - third dimension of the three dimensional array
3265: . mstart - first index you will use in first coordinate direction (often 0)
3266: . nstart - first index in the second coordinate direction (often 0)
3267: . pstart - first index in the third coordinate direction (often 0)
3268: - a - location of pointer to array obtained from VecGetArray3d()
3270: Level: developer
3272: Notes:
3273: For regular PETSc vectors this routine does not involve any copies. For
3274: any special vectors that do not store local vector data in a contiguous
3275: array, this routine will copy the data back into the underlying
3276: vector data structure from the array obtained with `VecGetArray()`.
3278: This routine actually zeros out the `a` pointer.
3280: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3281: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3282: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3283: @*/
3284: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3285: {
3286: void *dummy;
3288: PetscFunctionBegin;
3290: PetscAssertPointer(a, 8);
3292: dummy = (void *)(*a + mstart);
3293: PetscCall(PetscFree(dummy));
3294: PetscCall(VecRestoreArrayWrite(x, NULL));
3295: *a = NULL;
3296: PetscFunctionReturn(PETSC_SUCCESS);
3297: }
3299: /*@C
3300: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3301: processor's portion of the vector data. You MUST call `VecRestoreArray4d()`
3302: when you no longer need access to the array.
3304: Logically Collective
3306: Input Parameters:
3307: + x - the vector
3308: . m - first dimension of four dimensional array
3309: . n - second dimension of four dimensional array
3310: . p - third dimension of four dimensional array
3311: . q - fourth dimension of four dimensional array
3312: . mstart - first index you will use in first coordinate direction (often 0)
3313: . nstart - first index in the second coordinate direction (often 0)
3314: . pstart - first index in the third coordinate direction (often 0)
3315: - qstart - first index in the fourth coordinate direction (often 0)
3317: Output Parameter:
3318: . a - location to put pointer to the array
3320: Level: developer
3322: Notes:
3323: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3324: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3325: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3326: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3328: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3330: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3331: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3332: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3333: @*/
3334: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3335: {
3336: PetscInt i, N, j, k;
3337: PetscScalar *aa, ***b, **c;
3339: PetscFunctionBegin;
3341: PetscAssertPointer(a, 10);
3343: PetscCall(VecGetLocalSize(x, &N));
3344: 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);
3345: PetscCall(VecGetArray(x, &aa));
3347: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3348: b = (PetscScalar ***)((*a) + m);
3349: c = (PetscScalar **)(b + m * n);
3350: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3351: for (i = 0; i < m; i++)
3352: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3353: for (i = 0; i < m; i++)
3354: for (j = 0; j < n; j++)
3355: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3356: *a -= mstart;
3357: PetscFunctionReturn(PETSC_SUCCESS);
3358: }
3360: /*@C
3361: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3362: processor's portion of the vector data. You MUST call `VecRestoreArray4dWrite()`
3363: when you no longer need access to the array.
3365: Logically Collective
3367: Input Parameters:
3368: + x - the vector
3369: . m - first dimension of four dimensional array
3370: . n - second dimension of four dimensional array
3371: . p - third dimension of four dimensional array
3372: . q - fourth dimension of four dimensional array
3373: . mstart - first index you will use in first coordinate direction (often 0)
3374: . nstart - first index in the second coordinate direction (often 0)
3375: . pstart - first index in the third coordinate direction (often 0)
3376: - qstart - first index in the fourth coordinate direction (often 0)
3378: Output Parameter:
3379: . a - location to put pointer to the array
3381: Level: developer
3383: Notes:
3384: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3385: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3386: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3387: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3389: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3391: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3392: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3393: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3394: @*/
3395: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3396: {
3397: PetscInt i, N, j, k;
3398: PetscScalar *aa, ***b, **c;
3400: PetscFunctionBegin;
3402: PetscAssertPointer(a, 10);
3404: PetscCall(VecGetLocalSize(x, &N));
3405: 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);
3406: PetscCall(VecGetArrayWrite(x, &aa));
3408: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3409: b = (PetscScalar ***)((*a) + m);
3410: c = (PetscScalar **)(b + m * n);
3411: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3412: for (i = 0; i < m; i++)
3413: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3414: for (i = 0; i < m; i++)
3415: for (j = 0; j < n; j++)
3416: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3417: *a -= mstart;
3418: PetscFunctionReturn(PETSC_SUCCESS);
3419: }
3421: /*@C
3422: VecRestoreArray4d - Restores a vector after `VecGetArray4d()` has been called.
3424: Logically Collective
3426: Input Parameters:
3427: + x - the vector
3428: . m - first dimension of four dimensional array
3429: . n - second dimension of the four dimensional array
3430: . p - third dimension of the four dimensional array
3431: . q - fourth dimension of the four dimensional array
3432: . mstart - first index you will use in first coordinate direction (often 0)
3433: . nstart - first index in the second coordinate direction (often 0)
3434: . pstart - first index in the third coordinate direction (often 0)
3435: . qstart - first index in the fourth coordinate direction (often 0)
3436: - a - location of pointer to array obtained from VecGetArray4d()
3438: Level: developer
3440: Notes:
3441: For regular PETSc vectors this routine does not involve any copies. For
3442: any special vectors that do not store local vector data in a contiguous
3443: array, this routine will copy the data back into the underlying
3444: vector data structure from the array obtained with `VecGetArray()`.
3446: This routine actually zeros out the `a` pointer.
3448: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3449: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3450: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecGet`
3451: @*/
3452: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3453: {
3454: void *dummy;
3456: PetscFunctionBegin;
3458: PetscAssertPointer(a, 10);
3460: dummy = (void *)(*a + mstart);
3461: PetscCall(PetscFree(dummy));
3462: PetscCall(VecRestoreArray(x, NULL));
3463: *a = NULL;
3464: PetscFunctionReturn(PETSC_SUCCESS);
3465: }
3467: /*@C
3468: VecRestoreArray4dWrite - Restores a vector after `VecGetArray4dWrite()` has been called.
3470: Logically Collective
3472: Input Parameters:
3473: + x - the vector
3474: . m - first dimension of four dimensional array
3475: . n - second dimension of the four dimensional array
3476: . p - third dimension of the four dimensional array
3477: . q - fourth dimension of the four dimensional array
3478: . mstart - first index you will use in first coordinate direction (often 0)
3479: . nstart - first index in the second coordinate direction (often 0)
3480: . pstart - first index in the third coordinate direction (often 0)
3481: . qstart - first index in the fourth coordinate direction (often 0)
3482: - a - location of pointer to array obtained from `VecGetArray4d()`
3484: Level: developer
3486: Notes:
3487: For regular PETSc vectors this routine does not involve any copies. For
3488: any special vectors that do not store local vector data in a contiguous
3489: array, this routine will copy the data back into the underlying
3490: vector data structure from the array obtained with `VecGetArray()`.
3492: This routine actually zeros out the `a` pointer.
3494: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3495: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3496: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3497: @*/
3498: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3499: {
3500: void *dummy;
3502: PetscFunctionBegin;
3504: PetscAssertPointer(a, 10);
3506: dummy = (void *)(*a + mstart);
3507: PetscCall(PetscFree(dummy));
3508: PetscCall(VecRestoreArrayWrite(x, NULL));
3509: *a = NULL;
3510: PetscFunctionReturn(PETSC_SUCCESS);
3511: }
3513: /*@C
3514: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3515: processor's portion of the vector data. You MUST call `VecRestoreArray2dRead()`
3516: when you no longer need access to the array.
3518: Logically Collective
3520: Input Parameters:
3521: + x - the vector
3522: . m - first dimension of two dimensional array
3523: . n - second dimension of two dimensional array
3524: . mstart - first index you will use in first coordinate direction (often 0)
3525: - nstart - first index in the second coordinate direction (often 0)
3527: Output Parameter:
3528: . a - location to put pointer to the array
3530: Level: developer
3532: Notes:
3533: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
3534: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3535: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3536: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
3538: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3540: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3541: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3542: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3543: @*/
3544: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3545: {
3546: PetscInt i, N;
3547: const PetscScalar *aa;
3549: PetscFunctionBegin;
3551: PetscAssertPointer(a, 6);
3553: PetscCall(VecGetLocalSize(x, &N));
3554: 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);
3555: PetscCall(VecGetArrayRead(x, &aa));
3557: PetscCall(PetscMalloc1(m, a));
3558: for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3559: *a -= mstart;
3560: PetscFunctionReturn(PETSC_SUCCESS);
3561: }
3563: /*@C
3564: VecRestoreArray2dRead - Restores a vector after `VecGetArray2dRead()` has been called.
3566: Logically Collective
3568: Input Parameters:
3569: + x - the vector
3570: . m - first dimension of two dimensional array
3571: . n - second dimension of the two dimensional array
3572: . mstart - first index you will use in first coordinate direction (often 0)
3573: . nstart - first index in the second coordinate direction (often 0)
3574: - a - location of pointer to array obtained from VecGetArray2d()
3576: Level: developer
3578: Notes:
3579: For regular PETSc vectors this routine does not involve any copies. For
3580: any special vectors that do not store local vector data in a contiguous
3581: array, this routine will copy the data back into the underlying
3582: vector data structure from the array obtained with `VecGetArray()`.
3584: This routine actually zeros out the `a` pointer.
3586: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3587: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3588: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3589: @*/
3590: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3591: {
3592: void *dummy;
3594: PetscFunctionBegin;
3596: PetscAssertPointer(a, 6);
3598: dummy = (void *)(*a + mstart);
3599: PetscCall(PetscFree(dummy));
3600: PetscCall(VecRestoreArrayRead(x, NULL));
3601: *a = NULL;
3602: PetscFunctionReturn(PETSC_SUCCESS);
3603: }
3605: /*@C
3606: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3607: processor's portion of the vector data. You MUST call `VecRestoreArray1dRead()`
3608: when you no longer need access to the array.
3610: Logically Collective
3612: Input Parameters:
3613: + x - the vector
3614: . m - first dimension of two dimensional array
3615: - mstart - first index you will use in first coordinate direction (often 0)
3617: Output Parameter:
3618: . a - location to put pointer to the array
3620: Level: developer
3622: Notes:
3623: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3624: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3625: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3627: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3629: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3630: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3631: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3632: @*/
3633: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3634: {
3635: PetscInt N;
3637: PetscFunctionBegin;
3639: PetscAssertPointer(a, 4);
3641: PetscCall(VecGetLocalSize(x, &N));
3642: 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);
3643: PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3644: *a -= mstart;
3645: PetscFunctionReturn(PETSC_SUCCESS);
3646: }
3648: /*@C
3649: VecRestoreArray1dRead - Restores a vector after `VecGetArray1dRead()` has been called.
3651: Logically Collective
3653: Input Parameters:
3654: + x - the vector
3655: . m - first dimension of two dimensional array
3656: . mstart - first index you will use in first coordinate direction (often 0)
3657: - a - location of pointer to array obtained from `VecGetArray1dRead()`
3659: Level: developer
3661: Notes:
3662: For regular PETSc vectors this routine does not involve any copies. For
3663: any special vectors that do not store local vector data in a contiguous
3664: array, this routine will copy the data back into the underlying
3665: vector data structure from the array obtained with `VecGetArray1dRead()`.
3667: This routine actually zeros out the `a` pointer.
3669: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3670: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3671: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3672: @*/
3673: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3674: {
3675: PetscFunctionBegin;
3678: PetscCall(VecRestoreArrayRead(x, NULL));
3679: *a = NULL;
3680: PetscFunctionReturn(PETSC_SUCCESS);
3681: }
3683: /*@C
3684: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3685: processor's portion of the vector data. You MUST call `VecRestoreArray3dRead()`
3686: when you no longer need access to the array.
3688: Logically Collective
3690: Input Parameters:
3691: + x - the vector
3692: . m - first dimension of three dimensional array
3693: . n - second dimension of three dimensional array
3694: . p - third dimension of three dimensional array
3695: . mstart - first index you will use in first coordinate direction (often 0)
3696: . nstart - first index in the second coordinate direction (often 0)
3697: - pstart - first index in the third coordinate direction (often 0)
3699: Output Parameter:
3700: . a - location to put pointer to the array
3702: Level: developer
3704: Notes:
3705: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3706: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3707: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3708: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3dRead()`.
3710: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3712: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3713: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3714: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3715: @*/
3716: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3717: {
3718: PetscInt i, N, j;
3719: const PetscScalar *aa;
3720: PetscScalar **b;
3722: PetscFunctionBegin;
3724: PetscAssertPointer(a, 8);
3726: PetscCall(VecGetLocalSize(x, &N));
3727: 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);
3728: PetscCall(VecGetArrayRead(x, &aa));
3730: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3731: b = (PetscScalar **)((*a) + m);
3732: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3733: for (i = 0; i < m; i++)
3734: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset((PetscScalar *)aa, i * n * p + j * p - pstart);
3735: *a -= mstart;
3736: PetscFunctionReturn(PETSC_SUCCESS);
3737: }
3739: /*@C
3740: VecRestoreArray3dRead - Restores a vector after `VecGetArray3dRead()` has been called.
3742: Logically Collective
3744: Input Parameters:
3745: + x - the vector
3746: . m - first dimension of three dimensional array
3747: . n - second dimension of the three dimensional array
3748: . p - third dimension of the three dimensional array
3749: . mstart - first index you will use in first coordinate direction (often 0)
3750: . nstart - first index in the second coordinate direction (often 0)
3751: . pstart - first index in the third coordinate direction (often 0)
3752: - a - location of pointer to array obtained from `VecGetArray3dRead()`
3754: Level: developer
3756: Notes:
3757: For regular PETSc vectors this routine does not involve any copies. For
3758: any special vectors that do not store local vector data in a contiguous
3759: array, this routine will copy the data back into the underlying
3760: vector data structure from the array obtained with `VecGetArray()`.
3762: This routine actually zeros out the `a` pointer.
3764: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3765: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3766: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3767: @*/
3768: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3769: {
3770: void *dummy;
3772: PetscFunctionBegin;
3774: PetscAssertPointer(a, 8);
3776: dummy = (void *)(*a + mstart);
3777: PetscCall(PetscFree(dummy));
3778: PetscCall(VecRestoreArrayRead(x, NULL));
3779: *a = NULL;
3780: PetscFunctionReturn(PETSC_SUCCESS);
3781: }
3783: /*@C
3784: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3785: processor's portion of the vector data. You MUST call `VecRestoreArray4dRead()`
3786: when you no longer need access to the array.
3788: Logically Collective
3790: Input Parameters:
3791: + x - the vector
3792: . m - first dimension of four dimensional array
3793: . n - second dimension of four dimensional array
3794: . p - third dimension of four dimensional array
3795: . q - fourth dimension of four dimensional array
3796: . mstart - first index you will use in first coordinate direction (often 0)
3797: . nstart - first index in the second coordinate direction (often 0)
3798: . pstart - first index in the third coordinate direction (often 0)
3799: - qstart - first index in the fourth coordinate direction (often 0)
3801: Output Parameter:
3802: . a - location to put pointer to the array
3804: Level: beginner
3806: Notes:
3807: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3808: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3809: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3810: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3812: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3814: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3815: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3816: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3817: @*/
3818: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3819: {
3820: PetscInt i, N, j, k;
3821: const PetscScalar *aa;
3822: PetscScalar ***b, **c;
3824: PetscFunctionBegin;
3826: PetscAssertPointer(a, 10);
3828: PetscCall(VecGetLocalSize(x, &N));
3829: 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);
3830: PetscCall(VecGetArrayRead(x, &aa));
3832: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3833: b = (PetscScalar ***)((*a) + m);
3834: c = (PetscScalar **)(b + m * n);
3835: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3836: for (i = 0; i < m; i++)
3837: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3838: for (i = 0; i < m; i++)
3839: for (j = 0; j < n; j++)
3840: 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;
3841: *a -= mstart;
3842: PetscFunctionReturn(PETSC_SUCCESS);
3843: }
3845: /*@C
3846: VecRestoreArray4dRead - Restores a vector after `VecGetArray4d()` has been called.
3848: Logically Collective
3850: Input Parameters:
3851: + x - the vector
3852: . m - first dimension of four dimensional array
3853: . n - second dimension of the four dimensional array
3854: . p - third dimension of the four dimensional array
3855: . q - fourth dimension of the four dimensional array
3856: . mstart - first index you will use in first coordinate direction (often 0)
3857: . nstart - first index in the second coordinate direction (often 0)
3858: . pstart - first index in the third coordinate direction (often 0)
3859: . qstart - first index in the fourth coordinate direction (often 0)
3860: - a - location of pointer to array obtained from `VecGetArray4dRead()`
3862: Level: beginner
3864: Notes:
3865: For regular PETSc vectors this routine does not involve any copies. For
3866: any special vectors that do not store local vector data in a contiguous
3867: array, this routine will copy the data back into the underlying
3868: vector data structure from the array obtained with `VecGetArray()`.
3870: This routine actually zeros out the `a` pointer.
3872: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3873: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3874: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3875: @*/
3876: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3877: {
3878: void *dummy;
3880: PetscFunctionBegin;
3882: PetscAssertPointer(a, 10);
3884: dummy = (void *)(*a + mstart);
3885: PetscCall(PetscFree(dummy));
3886: PetscCall(VecRestoreArrayRead(x, NULL));
3887: *a = NULL;
3888: PetscFunctionReturn(PETSC_SUCCESS);
3889: }
3891: #if defined(PETSC_USE_DEBUG)
3893: /*@
3894: VecLockGet - Gets the current lock status of a vector
3896: Logically Collective
3898: Input Parameter:
3899: . x - the vector
3901: Output Parameter:
3902: . state - greater than zero indicates the vector is locked for read; less than zero indicates the vector is
3903: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3905: Level: advanced
3907: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3908: @*/
3909: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3910: {
3911: PetscFunctionBegin;
3913: *state = x->lock;
3914: PetscFunctionReturn(PETSC_SUCCESS);
3915: }
3917: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3918: {
3919: PetscFunctionBegin;
3921: PetscAssertPointer(file, 2);
3922: PetscAssertPointer(func, 3);
3923: PetscAssertPointer(line, 4);
3924: #if !PetscDefined(HAVE_THREADSAFETY)
3925: {
3926: const int index = x->lockstack.currentsize - 1;
3928: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupted vec lock stack, have negative index %d", index);
3929: *file = x->lockstack.file[index];
3930: *func = x->lockstack.function[index];
3931: *line = x->lockstack.line[index];
3932: }
3933: #else
3934: *file = NULL;
3935: *func = NULL;
3936: *line = 0;
3937: #endif
3938: PetscFunctionReturn(PETSC_SUCCESS);
3939: }
3941: /*@
3942: VecLockReadPush - Pushes a read-only lock on a vector to prevent it from being written to
3944: Logically Collective
3946: Input Parameter:
3947: . x - the vector
3949: Level: intermediate
3951: Notes:
3952: If this is set then calls to `VecGetArray()` or `VecSetValues()` or any other routines that change the vectors values will generate an error.
3954: The call can be nested, i.e., called multiple times on the same vector, but each `VecLockReadPush()` has to have one matching
3955: `VecLockReadPop()`, which removes the latest read-only lock.
3957: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3958: @*/
3959: PetscErrorCode VecLockReadPush(Vec x)
3960: {
3961: PetscFunctionBegin;
3963: 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");
3964: #if !PetscDefined(HAVE_THREADSAFETY)
3965: {
3966: const char *file, *func;
3967: int index, line;
3969: if ((index = petscstack.currentsize - 2) == -1) {
3970: // vec was locked "outside" of petsc, either in user-land or main. the error message will
3971: // now show this function as the culprit, but it will include the stacktrace
3972: file = "unknown user-file";
3973: func = "unknown_user_function";
3974: line = 0;
3975: } else {
3976: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected petscstack, have negative index %d", index);
3977: file = petscstack.file[index];
3978: func = petscstack.function[index];
3979: line = petscstack.line[index];
3980: }
3981: PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
3982: }
3983: #endif
3984: PetscFunctionReturn(PETSC_SUCCESS);
3985: }
3987: /*@
3988: VecLockReadPop - Pops a read-only lock from a vector
3990: Logically Collective
3992: Input Parameter:
3993: . x - the vector
3995: Level: intermediate
3997: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
3998: @*/
3999: PetscErrorCode VecLockReadPop(Vec x)
4000: {
4001: PetscFunctionBegin;
4003: PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
4004: #if !PetscDefined(HAVE_THREADSAFETY)
4005: {
4006: const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];
4008: PetscStackPop_Private(x->lockstack, previous);
4009: }
4010: #endif
4011: PetscFunctionReturn(PETSC_SUCCESS);
4012: }
4014: /*@C
4015: VecLockWriteSet - Lock or unlock a vector for exclusive read/write access
4017: Logically Collective
4019: Input Parameters:
4020: + x - the vector
4021: - flg - `PETSC_TRUE` to lock the vector for exclusive read/write access; `PETSC_FALSE` to unlock it.
4023: Level: intermediate
4025: Notes:
4026: The function is useful in split-phase computations, which usually have a begin phase and an end phase.
4027: One can call `VecLockWriteSet`(x,`PETSC_TRUE`) in the begin phase to lock a vector for exclusive
4028: access, and call `VecLockWriteSet`(x,`PETSC_FALSE`) in the end phase to unlock the vector from exclusive
4029: access. In this way, one is ensured no other operations can access the vector in between. The code may like
4031: .vb
4032: VecGetArray(x,&xdata); // begin phase
4033: VecLockWriteSet(v,PETSC_TRUE);
4035: Other operations, which can not access x anymore (they can access xdata, of course)
4037: VecRestoreArray(x,&vdata); // end phase
4038: VecLockWriteSet(v,PETSC_FALSE);
4039: .ve
4041: The call can not be nested on the same vector, in other words, one can not call `VecLockWriteSet`(x,`PETSC_TRUE`)
4042: again before calling `VecLockWriteSet`(v,`PETSC_FALSE`).
4044: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
4045: @*/
4046: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
4047: {
4048: PetscFunctionBegin;
4050: if (flg) {
4051: 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");
4052: 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");
4053: x->lock = -1;
4054: } else {
4055: 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");
4056: x->lock = 0;
4057: }
4058: PetscFunctionReturn(PETSC_SUCCESS);
4059: }
4061: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
4062: /*@
4063: VecLockPush - Pushes a read-only lock on a vector to prevent it from being written to
4065: Level: deprecated
4067: .seealso: [](ch_vectors), `Vec`, `VecLockReadPush()`
4068: @*/
4069: PetscErrorCode VecLockPush(Vec x)
4070: {
4071: PetscFunctionBegin;
4072: PetscCall(VecLockReadPush(x));
4073: PetscFunctionReturn(PETSC_SUCCESS);
4074: }
4076: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
4077: /*@
4078: VecLockPop - Pops a read-only lock from a vector
4080: Level: deprecated
4082: .seealso: [](ch_vectors), `Vec`, `VecLockReadPop()`
4083: @*/
4084: PetscErrorCode VecLockPop(Vec x)
4085: {
4086: PetscFunctionBegin;
4087: PetscCall(VecLockReadPop(x));
4088: PetscFunctionReturn(PETSC_SUCCESS);
4089: }
4091: #endif