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: Note:
92: For complex vectors, `VecDot()` computes
93: .vb
94: val = (x,y) = y^H x,
95: .ve
96: where $y^H$ denotes the conjugate transpose of `y`. Note that this corresponds to the usual "mathematicians" complex
97: inner product where the SECOND argument gets the complex conjugate. Since the `BLASdot()` complex conjugates the first
98: first argument we call the `BLASdot()` with the arguments reversed.
100: Use `VecTDot()` for the indefinite form
101: .vb
102: val = (x,y) = y^T x,
103: .ve
104: where $y^T$ denotes the transpose of `y`.
106: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDotRealPart()`
107: @*/
108: PetscErrorCode VecDot(Vec x, Vec y, PetscScalar *val)
109: {
110: PetscFunctionBegin;
113: PetscAssertPointer(val, 3);
116: PetscCheckSameTypeAndComm(x, 1, y, 2);
117: VecCheckSameSize(x, 1, y, 2);
118: VecCheckAssembled(x);
119: VecCheckAssembled(y);
121: PetscCall(VecLockReadPush(x));
122: PetscCall(VecLockReadPush(y));
123: PetscCall(PetscLogEventBegin(VEC_Dot, x, y, 0, 0));
124: PetscUseTypeMethod(x, dot, y, val);
125: PetscCall(PetscLogEventEnd(VEC_Dot, x, y, 0, 0));
126: PetscCall(VecLockReadPop(x));
127: PetscCall(VecLockReadPop(y));
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: /*@
132: VecDotRealPart - Computes the real part of the vector dot product.
134: Collective
136: Input Parameters:
137: + x - first vector
138: - y - second vector
140: Output Parameter:
141: . val - the real part of the dot product;
143: Level: intermediate
145: Notes for Users of Complex Numbers:
146: See `VecDot()` for more details on the definition of the dot product for complex numbers
148: For real numbers this returns the same value as `VecDot()`
150: 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
151: the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)
153: Developer Notes:
154: This is not currently optimized to compute only the real part of the dot product.
156: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDot()`, `VecDotNorm2()`
157: @*/
158: PetscErrorCode VecDotRealPart(Vec x, Vec y, PetscReal *val)
159: {
160: PetscScalar fdot;
162: PetscFunctionBegin;
163: PetscCall(VecDot(x, y, &fdot));
164: *val = PetscRealPart(fdot);
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: /*@
169: VecNorm - Computes the vector norm.
171: Collective
173: Input Parameters:
174: + x - the vector
175: - type - the type of the norm requested
177: Output Parameter:
178: . val - the norm
180: Level: intermediate
182: Notes:
183: See `NormType` for descriptions of each norm.
185: For complex numbers `NORM_1` will return the traditional 1 norm of the 2 norm of the complex
186: numbers; that is the 1 norm of the absolute values of the complex entries. In PETSc 3.6 and
187: earlier releases it returned the 1 norm of the 1 norm of the complex entries (what is
188: returned by the BLAS routine `asum()`). Both are valid norms but most people expect the former.
190: This routine stashes the computed norm value, repeated calls before the vector entries are
191: changed are then rapid since the precomputed value is immediately available. Certain vector
192: operations such as `VecSet()` store the norms so the value is immediately available and does
193: not need to be explicitly computed. `VecScale()` updates any stashed norm values, thus calls
194: after `VecScale()` do not need to explicitly recompute the norm.
196: .seealso: [](ch_vectors), `Vec`, `NormType`, `VecDot()`, `VecTDot()`, `VecDotBegin()`, `VecDotEnd()`, `VecNormAvailable()`,
197: `VecNormBegin()`, `VecNormEnd()`, `NormType()`
198: @*/
199: PetscErrorCode VecNorm(Vec x, NormType type, PetscReal *val)
200: {
201: PetscBool flg = PETSC_TRUE;
203: PetscFunctionBegin;
204: PetscCall(VecLockReadPush(x));
207: VecCheckAssembled(x);
209: PetscAssertPointer(val, 3);
211: PetscCall(VecNormAvailable(x, type, &flg, val));
212: // check that all MPI processes call this routine together and have same availability
213: if (PetscDefined(USE_DEBUG)) {
214: PetscMPIInt b0 = (PetscMPIInt)flg, b1[2], b2[2];
215: b1[0] = -b0;
216: b1[1] = b0;
217: PetscCallMPI(MPIU_Allreduce(b1, b2, 2, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
218: 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]);
219: if (flg) {
220: PetscReal b1[2], b2[2];
221: b1[0] = -(*val);
222: b1[1] = *val;
223: PetscCallMPI(MPIU_Allreduce(b1, b2, 2, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)x)));
224: PetscCheck((PetscIsNanReal(b2[0]) && PetscIsNanReal(b2[1])) || (-b2[0] == b2[1]), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Difference in cached %s norms: local %g", NormTypes[type], (double)*val);
225: }
226: }
227: if (!flg) {
228: PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
229: PetscUseTypeMethod(x, norm, type, val);
230: PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
232: if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
233: }
234: PetscCall(VecLockReadPop(x));
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: /*@
239: VecNormAvailable - Returns the vector norm if it is already known. That is, it has been previously computed and cached in the vector
241: Not Collective
243: Input Parameters:
244: + x - the vector
245: - 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
246: `NORM_1_AND_2`, which computes both norms and stores them
247: in a two element array.
249: Output Parameters:
250: + available - `PETSC_TRUE` if the val returned is valid
251: - val - the norm
253: Level: intermediate
255: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
256: `VecNormBegin()`, `VecNormEnd()`
257: @*/
258: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
259: {
260: PetscFunctionBegin;
263: PetscAssertPointer(available, 3);
264: PetscAssertPointer(val, 4);
266: if (type == NORM_1_AND_2) {
267: *available = PETSC_FALSE;
268: } else {
269: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
270: }
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: /*@
275: VecNormalize - Normalizes a vector by its 2-norm.
277: Collective
279: Input Parameter:
280: . x - the vector
282: Output Parameter:
283: . val - the vector norm before normalization. May be `NULL` if the value is not needed.
285: Level: intermediate
287: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
288: @*/
289: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
290: {
291: PetscReal norm;
293: PetscFunctionBegin;
296: PetscCall(VecSetErrorIfLocked(x, 1));
297: if (val) PetscAssertPointer(val, 2);
298: PetscCall(PetscLogEventBegin(VEC_Normalize, x, 0, 0, 0));
299: PetscCall(VecNorm(x, NORM_2, &norm));
300: if (norm == 0.0) PetscCall(PetscInfo(x, "Vector of zero norm can not be normalized; Returning only the zero norm\n"));
301: else if (PetscIsInfOrNanReal(norm)) PetscCall(PetscInfo(x, "Vector with infinity or NaN norm can not be normalized; Returning only the norm\n"));
302: else {
303: PetscScalar s = 1.0 / norm;
304: PetscCall(VecScale(x, s));
305: }
306: PetscCall(PetscLogEventEnd(VEC_Normalize, x, 0, 0, 0));
307: if (val) *val = norm;
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: /*@
312: VecMax - Determines the vector component with maximum real part and its location.
314: Collective
316: Input Parameter:
317: . x - the vector
319: Output Parameters:
320: + p - the index of `val` (pass `NULL` if you don't want this) in the vector
321: - val - the maximum component
323: Level: intermediate
325: Notes:
326: Returns the value `PETSC_MIN_REAL` and negative `p` if the vector is of length 0.
328: Returns the smallest index with the maximum value
330: Developer Note:
331: The Nag Fortran compiler does not like the symbol name VecMax
333: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
334: @*/
335: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
336: {
337: PetscFunctionBegin;
340: VecCheckAssembled(x);
341: if (p) PetscAssertPointer(p, 2);
342: PetscAssertPointer(val, 3);
343: PetscCall(VecLockReadPush(x));
344: PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
345: PetscUseTypeMethod(x, max, p, val);
346: PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
347: PetscCall(VecLockReadPop(x));
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: /*@
352: VecMin - Determines the vector component with minimum real part and its location.
354: Collective
356: Input Parameter:
357: . x - the vector
359: Output Parameters:
360: + p - the index of `val` (pass `NULL` if you don't want this location) in the vector
361: - val - the minimum component
363: Level: intermediate
365: Notes:
366: Returns the value `PETSC_MAX_REAL` and negative `p` if the vector is of length 0.
368: This returns the smallest index with the minimum value
370: Developer Note:
371: The Nag Fortran compiler does not like the symbol name VecMin
373: .seealso: [](ch_vectors), `Vec`, `VecMax()`
374: @*/
375: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
376: {
377: PetscFunctionBegin;
380: VecCheckAssembled(x);
381: if (p) PetscAssertPointer(p, 2);
382: PetscAssertPointer(val, 3);
383: PetscCall(VecLockReadPush(x));
384: PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
385: PetscUseTypeMethod(x, min, p, val);
386: PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
387: PetscCall(VecLockReadPop(x));
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: /*@
392: VecTDot - Computes an indefinite vector dot product. That is, this
393: routine does NOT use the complex conjugate.
395: Collective
397: Input Parameters:
398: + x - first vector
399: - y - second vector
401: Output Parameter:
402: . val - the dot product
404: Level: intermediate
406: Notes for Users of Complex Numbers:
407: For complex vectors, `VecTDot()` computes the indefinite form
408: .vb
409: val = (x,y) = y^T x,
410: .ve
411: where y^T denotes the transpose of y.
413: Use `VecDot()` for the inner product
414: .vb
415: val = (x,y) = y^H x,
416: .ve
417: where y^H denotes the conjugate transpose of y.
419: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
420: @*/
421: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
422: {
423: PetscFunctionBegin;
426: PetscAssertPointer(val, 3);
429: PetscCheckSameTypeAndComm(x, 1, y, 2);
430: VecCheckSameSize(x, 1, y, 2);
431: VecCheckAssembled(x);
432: VecCheckAssembled(y);
434: PetscCall(VecLockReadPush(x));
435: PetscCall(VecLockReadPush(y));
436: PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
437: PetscUseTypeMethod(x, tdot, y, val);
438: PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
439: PetscCall(VecLockReadPop(x));
440: PetscCall(VecLockReadPop(y));
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
444: PetscErrorCode VecScaleAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
445: {
446: PetscReal norms[4];
447: PetscBool flgs[4];
448: PetscScalar one = 1.0;
450: PetscFunctionBegin;
453: VecCheckAssembled(x);
454: PetscCall(VecSetErrorIfLocked(x, 1));
456: if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);
458: /* get current stashed norms */
459: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
461: PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
462: VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
463: PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));
465: PetscCall(PetscObjectStateIncrease((PetscObject)x));
466: /* put the scaled stashed norms back into the Vec */
467: for (PetscInt i = 0; i < 4; i++) {
468: PetscReal ar = PetscAbsScalar(alpha);
469: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
470: }
471: PetscFunctionReturn(PETSC_SUCCESS);
472: }
474: /*@
475: VecScale - Scales a vector.
477: Logically Collective
479: Input Parameters:
480: + x - the vector
481: - alpha - the scalar
483: Level: intermediate
485: Note:
486: For a vector with n components, `VecScale()` computes x[i] = alpha * x[i], for i=1,...,n.
488: .seealso: [](ch_vectors), `Vec`, `VecSet()`
489: @*/
490: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
491: {
492: PetscFunctionBegin;
493: PetscCall(VecScaleAsync_Private(x, alpha, NULL));
494: PetscFunctionReturn(PETSC_SUCCESS);
495: }
497: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
498: {
499: PetscFunctionBegin;
502: VecCheckAssembled(x);
504: PetscCall(VecSetErrorIfLocked(x, 1));
506: if (alpha == 0) {
507: PetscReal norm;
508: PetscBool set;
510: PetscCall(VecNormAvailable(x, NORM_2, &set, &norm));
511: if (set == PETSC_TRUE && norm == 0) PetscFunctionReturn(PETSC_SUCCESS);
512: }
513: PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
514: VecMethodDispatch(x, dctx, VecAsyncFnName(Set), set, (Vec, PetscScalar, PetscDeviceContext), alpha);
515: PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
516: PetscCall(PetscObjectStateIncrease((PetscObject)x));
518: /* norms can be simply set (if |alpha|*N not too large) */
519: {
520: PetscReal val = PetscAbsScalar(alpha);
521: const PetscInt N = x->map->N;
523: if (N == 0) {
524: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0));
525: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
526: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
527: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
528: } else if (val > PETSC_MAX_REAL / N) {
529: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
530: } else {
531: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
532: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
533: val *= PetscSqrtReal((PetscReal)N);
534: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
535: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
536: }
537: }
538: PetscFunctionReturn(PETSC_SUCCESS);
539: }
541: /*@
542: VecSet - Sets all components of a vector to a single scalar value.
544: Logically Collective
546: Input Parameters:
547: + x - the vector
548: - alpha - the scalar
550: Level: beginner
552: Notes:
553: For a vector of dimension n, `VecSet()` sets x[i] = alpha, for i=1,...,n,
554: so that all vector entries then equal the identical
555: scalar value, `alpha`. Use the more general routine
556: `VecSetValues()` to set different vector entries.
558: You CANNOT call this after you have called `VecSetValues()` but before you call
559: `VecAssemblyBegin()`
561: If `alpha` is zero and the norm of the vector is known to be zero then this skips the unneeded zeroing process
563: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
564: @*/
565: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
566: {
567: PetscFunctionBegin;
568: PetscCall(VecSetAsync_Private(x, alpha, NULL));
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: PetscErrorCode VecAXPYAsync_Private(Vec y, PetscScalar alpha, Vec x, PetscDeviceContext dctx)
573: {
574: PetscFunctionBegin;
579: PetscCheckSameTypeAndComm(x, 3, y, 1);
580: VecCheckSameSize(x, 3, y, 1);
581: VecCheckAssembled(x);
582: VecCheckAssembled(y);
584: if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
585: PetscCall(VecSetErrorIfLocked(y, 1));
586: if (x == y) {
587: PetscCall(VecScale(y, alpha + 1.0));
588: PetscFunctionReturn(PETSC_SUCCESS);
589: }
590: PetscCall(VecLockReadPush(x));
591: PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
592: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPY), axpy, (Vec, PetscScalar, Vec, PetscDeviceContext), alpha, x);
593: PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
594: PetscCall(VecLockReadPop(x));
595: PetscCall(PetscObjectStateIncrease((PetscObject)y));
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
598: /*@
599: VecAXPY - Computes `y = alpha x + y`.
601: Logically Collective
603: Input Parameters:
604: + alpha - the scalar
605: . x - vector scale by `alpha`
606: - y - vector accumulated into
608: Output Parameter:
609: . y - output vector
611: Level: intermediate
613: Notes:
614: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
615: .vb
616: VecAXPY(y,alpha,x) y = alpha x + y
617: VecAYPX(y,beta,x) y = x + beta y
618: VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
619: VecWAXPY(w,alpha,x,y) w = alpha x + y
620: VecAXPBYPCZ(z,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
621: VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
622: .ve
624: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
625: @*/
626: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
627: {
628: PetscFunctionBegin;
629: PetscCall(VecAXPYAsync_Private(y, alpha, x, NULL));
630: PetscFunctionReturn(PETSC_SUCCESS);
631: }
633: PetscErrorCode VecAYPXAsync_Private(Vec y, PetscScalar beta, Vec x, PetscDeviceContext dctx)
634: {
635: PetscFunctionBegin;
640: PetscCheckSameTypeAndComm(x, 3, y, 1);
641: VecCheckSameSize(x, 1, y, 3);
642: VecCheckAssembled(x);
643: VecCheckAssembled(y);
645: PetscCall(VecSetErrorIfLocked(y, 1));
646: if (x == y) {
647: PetscCall(VecScale(y, beta + 1.0));
648: PetscFunctionReturn(PETSC_SUCCESS);
649: }
650: PetscCall(VecLockReadPush(x));
651: if (beta == (PetscScalar)0.0) {
652: PetscCall(VecCopy(x, y));
653: } else {
654: PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
655: VecMethodDispatch(y, dctx, VecAsyncFnName(AYPX), aypx, (Vec, PetscScalar, Vec, PetscDeviceContext), beta, x);
656: PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
657: PetscCall(PetscObjectStateIncrease((PetscObject)y));
658: }
659: PetscCall(VecLockReadPop(x));
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
663: /*@
664: VecAYPX - Computes `y = x + beta y`.
666: Logically Collective
668: Input Parameters:
669: + beta - the scalar
670: . x - the unscaled vector
671: - y - the vector to be scaled
673: Output Parameter:
674: . y - output vector
676: Level: intermediate
678: Developer Notes:
679: The implementation is optimized for `beta` of -1.0, 0.0, and 1.0
681: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
682: @*/
683: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
684: {
685: PetscFunctionBegin;
686: PetscCall(VecAYPXAsync_Private(y, beta, x, NULL));
687: PetscFunctionReturn(PETSC_SUCCESS);
688: }
690: PetscErrorCode VecAXPBYAsync_Private(Vec y, PetscScalar alpha, PetscScalar beta, Vec x, PetscDeviceContext dctx)
691: {
692: PetscFunctionBegin;
697: PetscCheckSameTypeAndComm(x, 4, y, 1);
698: VecCheckSameSize(y, 1, x, 4);
699: VecCheckAssembled(x);
700: VecCheckAssembled(y);
703: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
704: if (x == y) {
705: PetscCall(VecScale(y, alpha + beta));
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
709: PetscCall(VecSetErrorIfLocked(y, 1));
710: PetscCall(VecLockReadPush(x));
711: PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
712: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPBY), axpby, (Vec, PetscScalar, PetscScalar, Vec, PetscDeviceContext), alpha, beta, x);
713: PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
714: PetscCall(PetscObjectStateIncrease((PetscObject)y));
715: PetscCall(VecLockReadPop(x));
716: PetscFunctionReturn(PETSC_SUCCESS);
717: }
719: /*@
720: VecAXPBY - Computes `y = alpha x + beta y`.
722: Logically Collective
724: Input Parameters:
725: + alpha - first scalar
726: . beta - second scalar
727: . x - the first scaled vector
728: - y - the second scaled vector
730: Output Parameter:
731: . y - output vector
733: Level: intermediate
735: Developer Notes:
736: The implementation is optimized for `alpha` and/or `beta` values of 0.0 and 1.0
738: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
739: @*/
740: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
741: {
742: PetscFunctionBegin;
743: PetscCall(VecAXPBYAsync_Private(y, alpha, beta, x, NULL));
744: PetscFunctionReturn(PETSC_SUCCESS);
745: }
747: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
748: {
749: PetscFunctionBegin;
756: PetscCheckSameTypeAndComm(x, 5, y, 6);
757: PetscCheckSameTypeAndComm(x, 5, z, 1);
758: VecCheckSameSize(x, 5, y, 6);
759: VecCheckSameSize(x, 5, z, 1);
760: PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
761: PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
762: VecCheckAssembled(x);
763: VecCheckAssembled(y);
764: VecCheckAssembled(z);
768: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
770: PetscCall(VecSetErrorIfLocked(z, 1));
771: PetscCall(VecLockReadPush(x));
772: PetscCall(VecLockReadPush(y));
773: PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
774: VecMethodDispatch(z, dctx, VecAsyncFnName(AXPBYPCZ), axpbypcz, (Vec, PetscScalar, PetscScalar, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, beta, gamma, x, y);
775: PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
776: PetscCall(PetscObjectStateIncrease((PetscObject)z));
777: PetscCall(VecLockReadPop(x));
778: PetscCall(VecLockReadPop(y));
779: PetscFunctionReturn(PETSC_SUCCESS);
780: }
781: /*@
782: VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`
784: Logically Collective
786: Input Parameters:
787: + alpha - first scalar
788: . beta - second scalar
789: . gamma - third scalar
790: . x - first vector
791: . y - second vector
792: - z - third vector
794: Output Parameter:
795: . z - output vector
797: Level: intermediate
799: Note:
800: `x`, `y` and `z` must be different vectors
802: Developer Notes:
803: The implementation is optimized for `alpha` of 1.0 and `gamma` of 1.0 or 0.0
805: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
806: @*/
807: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
808: {
809: PetscFunctionBegin;
810: PetscCall(VecAXPBYPCZAsync_Private(z, alpha, beta, gamma, x, y, NULL));
811: PetscFunctionReturn(PETSC_SUCCESS);
812: }
814: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
815: {
816: PetscFunctionBegin;
823: PetscCheckSameTypeAndComm(x, 3, y, 4);
824: PetscCheckSameTypeAndComm(y, 4, w, 1);
825: VecCheckSameSize(x, 3, y, 4);
826: VecCheckSameSize(x, 3, w, 1);
827: PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
828: PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
829: VecCheckAssembled(x);
830: VecCheckAssembled(y);
832: PetscCall(VecSetErrorIfLocked(w, 1));
834: PetscCall(VecLockReadPush(x));
835: PetscCall(VecLockReadPush(y));
836: if (alpha == (PetscScalar)0.0) {
837: PetscCall(VecCopyAsync_Private(y, w, dctx));
838: } else {
839: PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
840: VecMethodDispatch(w, dctx, VecAsyncFnName(WAXPY), waxpy, (Vec, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, x, y);
841: PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
842: PetscCall(PetscObjectStateIncrease((PetscObject)w));
843: }
844: PetscCall(VecLockReadPop(x));
845: PetscCall(VecLockReadPop(y));
846: PetscFunctionReturn(PETSC_SUCCESS);
847: }
849: /*@
850: VecWAXPY - Computes `w = alpha x + y`.
852: Logically Collective
854: Input Parameters:
855: + alpha - the scalar
856: . x - first vector, multiplied by `alpha`
857: - y - second vector
859: Output Parameter:
860: . w - the result
862: Level: intermediate
864: Note:
865: `w` cannot be either `x` or `y`, but `x` and `y` can be the same
867: Developer Notes:
868: The implementation is optimized for alpha of -1.0, 0.0, and 1.0
870: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
871: @*/
872: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
873: {
874: PetscFunctionBegin;
875: PetscCall(VecWAXPYAsync_Private(w, alpha, x, y, NULL));
876: PetscFunctionReturn(PETSC_SUCCESS);
877: }
879: /*@
880: VecSetValues - Inserts or adds values into certain locations of a vector.
882: Not Collective
884: Input Parameters:
885: + x - vector to insert in
886: . ni - number of elements to add
887: . ix - indices where to add
888: . y - array of values. Pass `NULL` to set all zeroes.
889: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries
891: Level: beginner
893: Notes:
894: .vb
895: `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
896: .ve
898: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
899: options cannot be mixed without intervening calls to the assembly
900: routines.
902: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
903: MUST be called after all calls to `VecSetValues()` have been completed.
905: VecSetValues() uses 0-based indices in Fortran as well as in C.
907: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
908: negative indices may be passed in ix. These rows are
909: simply ignored. This allows easily inserting element load matrices
910: with homogeneous Dirichlet boundary conditions that you don't want represented
911: in the vector.
913: Fortran Note:
914: If any of `ix` and `y` are scalars pass them using, for example,
915: .vb
916: call VecSetValues(mat, one, [ix], [y], INSERT_VALUES, ierr)
917: .ve
919: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
920: `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`,
921: `VecOption`, `VecSetOption()`
922: @*/
923: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
924: {
925: PetscFunctionBeginHot;
927: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
928: PetscAssertPointer(ix, 3);
929: if (y) PetscAssertPointer(y, 4);
932: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
933: PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
934: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
935: PetscCall(PetscObjectStateIncrease((PetscObject)x));
936: PetscFunctionReturn(PETSC_SUCCESS);
937: }
939: /*@
940: VecGetValues - Gets values from certain locations of a vector. Currently
941: can only get values on the same processor on which they are owned
943: Not Collective
945: Input Parameters:
946: + x - vector to get values from
947: . ni - number of elements to get
948: - ix - indices where to get them from (in global 1d numbering)
950: Output Parameter:
951: . y - array of values, must be passed in with a length of `ni`
953: Level: beginner
955: Notes:
956: The user provides the allocated array y; it is NOT allocated in this routine
958: `VecGetValues()` gets y[i] = x[ix[i]], for i=0,...,ni-1.
960: `VecAssemblyBegin()` and `VecAssemblyEnd()` MUST be called before calling this if `VecSetValues()` or related routine has been called
962: VecGetValues() uses 0-based indices in Fortran as well as in C.
964: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
965: negative indices may be passed in ix. These rows are
966: simply ignored.
968: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
969: @*/
970: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
971: {
972: PetscFunctionBegin;
974: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
975: PetscAssertPointer(ix, 3);
976: PetscAssertPointer(y, 4);
978: VecCheckAssembled(x);
979: PetscUseTypeMethod(x, getvalues, ni, ix, y);
980: PetscFunctionReturn(PETSC_SUCCESS);
981: }
983: /*@
984: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
986: Not Collective
988: Input Parameters:
989: + x - vector to insert in
990: . ni - number of blocks to add
991: . ix - indices where to add in block count, rather than element count
992: . y - array of values. Pass `NULL` to set all zeroes.
993: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries
995: Level: intermediate
997: Notes:
998: `VecSetValuesBlocked()` sets x[bs*ix[i]+j] = y[bs*i+j],
999: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
1001: Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
1002: options cannot be mixed without intervening calls to the assembly
1003: routines.
1005: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1006: MUST be called after all calls to `VecSetValuesBlocked()` have been completed.
1008: `VecSetValuesBlocked()` uses 0-based indices in Fortran as well as in C.
1010: Negative indices may be passed in ix, these rows are
1011: simply ignored. This allows easily inserting element load matrices
1012: with homogeneous Dirichlet boundary conditions that you don't want represented
1013: in the vector.
1015: Fortran Note:
1016: If any of `ix` and `y` are scalars pass them using, for example,
1017: .vb
1018: call VecSetValuesBlocked(mat, one, [ix], [y], INSERT_VALUES, ierr)
1019: .ve
1021: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
1022: `VecSetValues()`
1023: @*/
1024: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1025: {
1026: PetscFunctionBeginHot;
1028: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1029: PetscAssertPointer(ix, 3);
1030: if (y) PetscAssertPointer(y, 4);
1033: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1034: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1035: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1036: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1037: PetscFunctionReturn(PETSC_SUCCESS);
1038: }
1040: /*@
1041: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1042: using a local ordering of the nodes.
1044: Not Collective
1046: Input Parameters:
1047: + x - vector to insert in
1048: . ni - number of elements to add
1049: . ix - indices where to add
1050: . y - array of values. Pass `NULL` to set all zeroes.
1051: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1053: Level: intermediate
1055: Notes:
1056: `VecSetValuesLocal()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
1058: Calls to `VecSetValuesLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1059: options cannot be mixed without intervening calls to the assembly
1060: routines.
1062: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1063: MUST be called after all calls to `VecSetValuesLocal()` have been completed.
1065: `VecSetValuesLocal()` uses 0-based indices in Fortran as well as in C.
1067: Fortran Note:
1068: If any of `ix` and `y` are scalars pass them using, for example,
1069: .vb
1070: call VecSetValuesLocal(mat, one, [ix], [y], INSERT_VALUES, ierr)
1071: .ve
1073: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1074: `VecSetValuesBlockedLocal()`
1075: @*/
1076: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1077: {
1078: PetscInt lixp[128], *lix = lixp;
1080: PetscFunctionBeginHot;
1082: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1083: PetscAssertPointer(ix, 3);
1084: if (y) PetscAssertPointer(y, 4);
1087: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1088: if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
1089: if (x->map->mapping) {
1090: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1091: PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1092: PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1093: if (ni > 128) PetscCall(PetscFree(lix));
1094: } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1095: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1096: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1097: PetscFunctionReturn(PETSC_SUCCESS);
1098: }
1100: /*@
1101: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1102: using a local ordering of the nodes.
1104: Not Collective
1106: Input Parameters:
1107: + x - vector to insert in
1108: . ni - number of blocks to add
1109: . ix - indices where to add in block count, not element count
1110: . y - array of values. Pass `NULL` to set all zeroes.
1111: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1113: Level: intermediate
1115: Notes:
1116: `VecSetValuesBlockedLocal()` sets x[bs*ix[i]+j] = y[bs*i+j],
1117: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with `VecSetBlockSize()`.
1119: Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1120: options cannot be mixed without intervening calls to the assembly
1121: routines.
1123: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1124: MUST be called after all calls to `VecSetValuesBlockedLocal()` have been completed.
1126: `VecSetValuesBlockedLocal()` uses 0-based indices in Fortran as well as in C.
1128: Fortran Note:
1129: If any of `ix` and `y` are scalars pass them using, for example,
1130: .vb
1131: call VecSetValuesBlockedLocal(mat, one, [ix], [y], INSERT_VALUES, ierr)
1132: .ve
1134: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1135: `VecSetLocalToGlobalMapping()`
1136: @*/
1137: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1138: {
1139: PetscInt lixp[128], *lix = lixp;
1141: PetscFunctionBeginHot;
1143: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1144: PetscAssertPointer(ix, 3);
1145: if (y) PetscAssertPointer(y, 4);
1147: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1148: if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
1149: if (x->map->mapping) {
1150: if (ni > (PetscInt)PETSC_STATIC_ARRAY_LENGTH(lixp)) PetscCall(PetscMalloc1(ni, &lix));
1151: PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1152: PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1153: if (ni > (PetscInt)PETSC_STATIC_ARRAY_LENGTH(lixp)) PetscCall(PetscFree(lix));
1154: } else {
1155: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1156: }
1157: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1158: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1159: PetscFunctionReturn(PETSC_SUCCESS);
1160: }
1162: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1163: {
1164: PetscFunctionBegin;
1167: VecCheckAssembled(x);
1169: if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1170: PetscAssertPointer(y, 3);
1171: for (PetscInt i = 0; i < nv; ++i) {
1174: PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1175: VecCheckSameSize(x, 1, y[i], 3);
1176: VecCheckAssembled(y[i]);
1177: PetscCall(VecLockReadPush(y[i]));
1178: }
1179: PetscAssertPointer(result, 4);
1182: PetscCall(VecLockReadPush(x));
1183: PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1184: PetscCall((*mxdot)(x, nv, y, result));
1185: PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1186: PetscCall(VecLockReadPop(x));
1187: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1188: PetscFunctionReturn(PETSC_SUCCESS);
1189: }
1191: /*@
1192: VecMTDot - Computes indefinite vector multiple dot products.
1193: That is, it does NOT use the complex conjugate.
1195: Collective
1197: Input Parameters:
1198: + x - one vector
1199: . nv - number of vectors
1200: - y - array of vectors. Note that vectors are pointers
1202: Output Parameter:
1203: . val - array of the dot products
1205: Level: intermediate
1207: Notes for Users of Complex Numbers:
1208: For complex vectors, `VecMTDot()` computes the indefinite form
1209: .vb
1210: val = (x,y) = y^T x,
1211: .ve
1212: where y^T denotes the transpose of y.
1214: Use `VecMDot()` for the inner product
1215: .vb
1216: val = (x,y) = y^H x,
1217: .ve
1218: where y^H denotes the conjugate transpose of y.
1220: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1221: @*/
1222: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1223: {
1224: PetscFunctionBegin;
1226: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1227: PetscFunctionReturn(PETSC_SUCCESS);
1228: }
1230: /*@
1231: VecMDot - Computes multiple vector dot products.
1233: Collective
1235: Input Parameters:
1236: + x - one vector
1237: . nv - number of vectors
1238: - y - array of vectors.
1240: Output Parameter:
1241: . val - array of the dot products (does not allocate the array)
1243: Level: intermediate
1245: Notes for Users of Complex Numbers:
1246: For complex vectors, `VecMDot()` computes
1247: .vb
1248: val = (x,y) = y^H x,
1249: .ve
1250: where y^H denotes the conjugate transpose of y.
1252: Use `VecMTDot()` for the indefinite form
1253: .vb
1254: val = (x,y) = y^T x,
1255: .ve
1256: where y^T denotes the transpose of y.
1258: Note:
1259: The implementation may use BLAS 2 operations when the vectors `y` have been obtained with `VecDuplicateVecs()`
1261: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`, `VecDuplicateVecs()`
1262: @*/
1263: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1264: {
1265: PetscFunctionBegin;
1267: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1268: PetscFunctionReturn(PETSC_SUCCESS);
1269: }
1271: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
1272: {
1273: PetscFunctionBegin;
1275: VecCheckAssembled(y);
1277: PetscCall(VecSetErrorIfLocked(y, 1));
1278: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1279: if (nv) {
1280: PetscInt zeros = 0;
1282: PetscAssertPointer(alpha, 3);
1283: PetscAssertPointer(x, 4);
1284: for (PetscInt i = 0; i < nv; ++i) {
1288: PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1289: VecCheckSameSize(y, 1, x[i], 4);
1290: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1291: VecCheckAssembled(x[i]);
1292: PetscCall(VecLockReadPush(x[i]));
1293: zeros += alpha[i] == (PetscScalar)0.0;
1294: }
1296: if (zeros < nv) {
1297: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1298: VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1299: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1300: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1301: }
1303: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1304: }
1305: PetscFunctionReturn(PETSC_SUCCESS);
1306: }
1308: /*@
1309: VecMAXPY - Computes `y = y + sum alpha[i] x[i]`
1311: Logically Collective
1313: Input Parameters:
1314: + nv - number of scalars and `x` vectors
1315: . alpha - array of scalars
1316: . y - one vector
1317: - x - array of vectors
1319: Level: intermediate
1321: Notes:
1322: `y` cannot be any of the `x` vectors
1324: The implementation may use BLAS 2 operations when the vectors `y` have been obtained with `VecDuplicateVecs()`
1326: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`, `VecDuplicateVecs()`
1327: @*/
1328: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1329: {
1330: PetscFunctionBegin;
1331: PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1332: PetscFunctionReturn(PETSC_SUCCESS);
1333: }
1335: /*@
1336: VecMAXPBY - Computes `y = beta y + sum alpha[i] x[i]`
1338: Logically Collective
1340: Input Parameters:
1341: + nv - number of scalars and `x` vectors
1342: . alpha - array of scalars
1343: . beta - scalar
1344: . y - one vector
1345: - x - array of vectors
1347: Level: intermediate
1349: Note:
1350: `y` cannot be any of the `x` vectors.
1352: Developer Notes:
1353: This is a convenience routine, but implementations might be able to optimize it, for example, when `beta` is zero.
1355: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1356: @*/
1357: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
1358: {
1359: PetscFunctionBegin;
1361: VecCheckAssembled(y);
1363: PetscCall(VecSetErrorIfLocked(y, 1));
1364: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1367: if (y->ops->maxpby) {
1368: PetscInt zeros = 0;
1370: if (nv) {
1371: PetscAssertPointer(alpha, 3);
1372: PetscAssertPointer(x, 5);
1373: }
1375: for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
1379: PetscCheckSameTypeAndComm(y, 1, x[i], 5);
1380: VecCheckSameSize(y, 1, x[i], 5);
1381: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1382: VecCheckAssembled(x[i]);
1383: PetscCall(VecLockReadPush(x[i]));
1384: zeros += alpha[i] == (PetscScalar)0.0;
1385: }
1387: if (zeros < nv) { // has nonzero alpha
1388: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1389: PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1390: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1391: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1392: } else {
1393: PetscCall(VecScale(y, beta));
1394: }
1396: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1397: } else { // no maxpby
1398: if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1399: else PetscCall(VecScale(y, beta));
1400: PetscCall(VecMAXPY(y, nv, alpha, x));
1401: }
1402: PetscFunctionReturn(PETSC_SUCCESS);
1403: }
1405: /*@
1406: VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
1407: in the order they appear in the array. The concatenated vector resides on the same
1408: communicator and is the same type as the source vectors.
1410: Collective
1412: Input Parameters:
1413: + nx - number of vectors to be concatenated
1414: - X - array containing the vectors to be concatenated in the order of concatenation
1416: Output Parameters:
1417: + Y - concatenated vector
1418: - x_is - array of index sets corresponding to the concatenated components of `Y` (pass `NULL` if not needed)
1420: Level: advanced
1422: Notes:
1423: Concatenation is similar to the functionality of a `VECNEST` object; they both represent combination of
1424: different vector spaces. However, concatenated vectors do not store any information about their
1425: sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1426: manipulation of data in the concatenated vector that corresponds to the original components at creation.
1428: This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1429: has to operate on combined vector spaces and cannot utilize `VECNEST` objects due to incompatibility with
1430: bound projections.
1432: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1433: @*/
1434: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1435: {
1436: MPI_Comm comm;
1437: VecType vec_type;
1438: Vec Ytmp, Xtmp;
1439: IS *is_tmp;
1440: PetscInt i, shift = 0, Xnl, Xng, Xbegin;
1442: PetscFunctionBegin;
1446: PetscAssertPointer(Y, 3);
1448: if ((*X)->ops->concatenate) {
1449: /* use the dedicated concatenation function if available */
1450: PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1451: } else {
1452: /* loop over vectors and start creating IS */
1453: comm = PetscObjectComm((PetscObject)*X);
1454: PetscCall(VecGetType(*X, &vec_type));
1455: PetscCall(PetscMalloc1(nx, &is_tmp));
1456: for (i = 0; i < nx; i++) {
1457: PetscCall(VecGetSize(X[i], &Xng));
1458: PetscCall(VecGetLocalSize(X[i], &Xnl));
1459: PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1460: PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1461: shift += Xng;
1462: }
1463: /* create the concatenated vector */
1464: PetscCall(VecCreate(comm, &Ytmp));
1465: PetscCall(VecSetType(Ytmp, vec_type));
1466: PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1467: PetscCall(VecSetUp(Ytmp));
1468: /* copy data from X array to Y and return */
1469: for (i = 0; i < nx; i++) {
1470: PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1471: PetscCall(VecCopy(X[i], Xtmp));
1472: PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1473: }
1474: *Y = Ytmp;
1475: if (x_is) {
1476: *x_is = is_tmp;
1477: } else {
1478: for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1479: PetscCall(PetscFree(is_tmp));
1480: }
1481: }
1482: PetscFunctionReturn(PETSC_SUCCESS);
1483: }
1485: /* A helper function for VecGetSubVector to check if we can implement it with no-copy (i.e. the subvector shares
1486: memory with the original vector), and the block size of the subvector.
1488: Input Parameters:
1489: + X - the original vector
1490: - is - the index set of the subvector
1492: Output Parameters:
1493: + contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
1494: . start - start of contiguous block, as an offset from the start of the ownership range of the original vector
1495: - blocksize - the block size of the subvector
1497: */
1498: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1499: {
1500: PetscInt gstart, gend, lstart;
1501: PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1502: PetscInt n, N, ibs, vbs, bs = 1;
1504: PetscFunctionBegin;
1505: PetscCall(ISGetLocalSize(is, &n));
1506: PetscCall(ISGetSize(is, &N));
1507: PetscCall(ISGetBlockSize(is, &ibs));
1508: PetscCall(VecGetBlockSize(X, &vbs));
1509: PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1510: PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1511: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1512: if (ibs > 1) {
1513: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1514: bs = ibs;
1515: } else {
1516: if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1517: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1518: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1519: }
1521: *contig = red[0];
1522: *start = lstart;
1523: *blocksize = bs;
1524: PetscFunctionReturn(PETSC_SUCCESS);
1525: }
1527: /* A helper function for VecGetSubVector, to be used when we have to build a standalone subvector through VecScatter
1529: Input Parameters:
1530: + X - the original vector
1531: . is - the index set of the subvector
1532: - bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()
1534: Output Parameter:
1535: . Z - the subvector, which will compose the VecScatter context on output
1536: */
1537: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1538: {
1539: PetscInt n, N;
1540: VecScatter vscat;
1541: Vec Y;
1543: PetscFunctionBegin;
1544: PetscCall(ISGetLocalSize(is, &n));
1545: PetscCall(ISGetSize(is, &N));
1546: PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1547: PetscCall(VecSetSizes(Y, n, N));
1548: PetscCall(VecSetBlockSize(Y, bs));
1549: PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1550: PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1551: PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1552: PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1553: PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1554: PetscCall(VecScatterDestroy(&vscat));
1555: *Z = Y;
1556: PetscFunctionReturn(PETSC_SUCCESS);
1557: }
1559: /*@
1560: VecGetSubVector - Gets a vector representing part of another vector
1562: Collective
1564: Input Parameters:
1565: + X - vector from which to extract a subvector
1566: - is - index set representing portion of `X` to extract
1568: Output Parameter:
1569: . Y - subvector corresponding to `is`
1571: Level: advanced
1573: Notes:
1574: The subvector `Y` should be returned with `VecRestoreSubVector()`.
1575: `X` and `is` must be defined on the same communicator
1577: Changes to the subvector will be reflected in the `X` vector on the call to `VecRestoreSubVector()`.
1579: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1580: modifying the subvector. Other non-overlapping subvectors can still be obtained from `X` using this function.
1582: 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`.
1584: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1585: @*/
1586: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1587: {
1588: Vec Z;
1590: PetscFunctionBegin;
1593: PetscCheckSameComm(X, 1, is, 2);
1594: PetscAssertPointer(Y, 3);
1595: if (X->ops->getsubvector) {
1596: PetscUseTypeMethod(X, getsubvector, is, &Z);
1597: } else { /* Default implementation currently does no caching */
1598: PetscBool contig;
1599: PetscInt n, N, start, bs;
1601: PetscCall(ISGetLocalSize(is, &n));
1602: PetscCall(ISGetSize(is, &N));
1603: PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1604: if (contig) { /* We can do a no-copy implementation */
1605: const PetscScalar *x;
1606: PetscInt state = 0;
1607: PetscBool isstd, iscuda, iship;
1609: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1610: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1611: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1612: if (iscuda) {
1613: #if defined(PETSC_HAVE_CUDA)
1614: const PetscScalar *x_d;
1615: PetscMPIInt size;
1616: PetscOffloadMask flg;
1618: PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1619: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1620: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1621: if (x) x += start;
1622: if (x_d) x_d += start;
1623: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1624: if (size == 1) {
1625: PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1626: } else {
1627: PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1628: }
1629: Z->offloadmask = flg;
1630: #endif
1631: } else if (iship) {
1632: #if defined(PETSC_HAVE_HIP)
1633: const PetscScalar *x_d;
1634: PetscMPIInt size;
1635: PetscOffloadMask flg;
1637: PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1638: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1639: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1640: if (x) x += start;
1641: if (x_d) x_d += start;
1642: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1643: if (size == 1) {
1644: PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1645: } else {
1646: PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1647: }
1648: Z->offloadmask = flg;
1649: #endif
1650: } else if (isstd) {
1651: PetscMPIInt size;
1653: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1654: PetscCall(VecGetArrayRead(X, &x));
1655: if (x) x += start;
1656: if (size == 1) {
1657: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1658: } else {
1659: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1660: }
1661: PetscCall(VecRestoreArrayRead(X, &x));
1662: } else { /* default implementation: use place array */
1663: PetscCall(VecGetArrayRead(X, &x));
1664: PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1665: PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1666: PetscCall(VecSetSizes(Z, n, N));
1667: PetscCall(VecSetBlockSize(Z, bs));
1668: PetscCall(VecPlaceArray(Z, PetscSafePointerPlusOffset(x, start)));
1669: PetscCall(VecRestoreArrayRead(X, &x));
1670: }
1672: /* this is relevant only in debug mode */
1673: PetscCall(VecLockGet(X, &state));
1674: if (state) PetscCall(VecLockReadPush(Z));
1675: Z->ops->placearray = NULL;
1676: Z->ops->replacearray = NULL;
1677: } else { /* Have to create a scatter and do a copy */
1678: PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1679: }
1680: }
1681: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1682: if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1683: PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1684: *Y = Z;
1685: PetscFunctionReturn(PETSC_SUCCESS);
1686: }
1688: /*@
1689: VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`
1691: Collective
1693: Input Parameters:
1694: + X - vector from which subvector was obtained
1695: . is - index set representing the subset of `X`
1696: - Y - subvector being restored
1698: Level: advanced
1700: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1701: @*/
1702: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1703: {
1704: PETSC_UNUSED PetscObjectState dummystate = 0;
1705: PetscBool unchanged;
1707: PetscFunctionBegin;
1710: PetscCheckSameComm(X, 1, is, 2);
1711: PetscAssertPointer(Y, 3);
1714: if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1715: else {
1716: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1717: if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1718: VecScatter scatter;
1719: PetscInt state;
1721: PetscCall(VecLockGet(X, &state));
1722: PetscCheck(state == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec X is locked for read-only or read/write access");
1724: PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1725: if (scatter) {
1726: PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1727: PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1728: } else {
1729: PetscBool iscuda, iship;
1730: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1731: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1733: if (iscuda) {
1734: #if defined(PETSC_HAVE_CUDA)
1735: PetscOffloadMask ymask = (*Y)->offloadmask;
1737: /* The offloadmask of X dictates where to move memory
1738: If X GPU data is valid, then move Y data on GPU if needed
1739: Otherwise, move back to the CPU */
1740: switch (X->offloadmask) {
1741: case PETSC_OFFLOAD_BOTH:
1742: if (ymask == PETSC_OFFLOAD_CPU) {
1743: PetscCall(VecCUDAResetArray(*Y));
1744: } else if (ymask == PETSC_OFFLOAD_GPU) {
1745: X->offloadmask = PETSC_OFFLOAD_GPU;
1746: }
1747: break;
1748: case PETSC_OFFLOAD_GPU:
1749: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1750: break;
1751: case PETSC_OFFLOAD_CPU:
1752: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1753: break;
1754: case PETSC_OFFLOAD_UNALLOCATED:
1755: case PETSC_OFFLOAD_KOKKOS:
1756: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1757: }
1758: #endif
1759: } else if (iship) {
1760: #if defined(PETSC_HAVE_HIP)
1761: PetscOffloadMask ymask = (*Y)->offloadmask;
1763: /* The offloadmask of X dictates where to move memory
1764: If X GPU data is valid, then move Y data on GPU if needed
1765: Otherwise, move back to the CPU */
1766: switch (X->offloadmask) {
1767: case PETSC_OFFLOAD_BOTH:
1768: if (ymask == PETSC_OFFLOAD_CPU) {
1769: PetscCall(VecHIPResetArray(*Y));
1770: } else if (ymask == PETSC_OFFLOAD_GPU) {
1771: X->offloadmask = PETSC_OFFLOAD_GPU;
1772: }
1773: break;
1774: case PETSC_OFFLOAD_GPU:
1775: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1776: break;
1777: case PETSC_OFFLOAD_CPU:
1778: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1779: break;
1780: case PETSC_OFFLOAD_UNALLOCATED:
1781: case PETSC_OFFLOAD_KOKKOS:
1782: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1783: }
1784: #endif
1785: } else {
1786: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1787: PetscCall(VecResetArray(*Y));
1788: }
1789: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1790: }
1791: }
1792: }
1793: PetscCall(VecDestroy(Y));
1794: PetscFunctionReturn(PETSC_SUCCESS);
1795: }
1797: /*@
1798: VecCreateLocalVector - Creates a vector object suitable for use with `VecGetLocalVector()` and friends. You must call `VecDestroy()` when the
1799: vector is no longer needed.
1801: Not Collective.
1803: Input Parameter:
1804: . v - The vector for which the local vector is desired.
1806: Output Parameter:
1807: . w - Upon exit this contains the local vector.
1809: Level: beginner
1811: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1812: @*/
1813: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1814: {
1815: VecType roottype;
1816: PetscInt n;
1818: PetscFunctionBegin;
1820: PetscAssertPointer(w, 2);
1821: if (v->ops->createlocalvector) {
1822: PetscUseTypeMethod(v, createlocalvector, w);
1823: PetscFunctionReturn(PETSC_SUCCESS);
1824: }
1825: PetscCall(VecGetRootType_Private(v, &roottype));
1826: PetscCall(VecCreate(PETSC_COMM_SELF, w));
1827: PetscCall(VecGetLocalSize(v, &n));
1828: PetscCall(VecSetSizes(*w, n, n));
1829: PetscCall(VecGetBlockSize(v, &n));
1830: PetscCall(VecSetBlockSize(*w, n));
1831: PetscCall(VecSetType(*w, roottype));
1832: PetscFunctionReturn(PETSC_SUCCESS);
1833: }
1835: /*@
1836: VecGetLocalVectorRead - Maps the local portion of a vector into a
1837: vector.
1839: Not Collective.
1841: Input Parameter:
1842: . v - The vector for which the local vector is desired.
1844: Output Parameter:
1845: . w - Upon exit this contains the local vector.
1847: Level: beginner
1849: Notes:
1850: You must call `VecRestoreLocalVectorRead()` when the local
1851: vector is no longer needed.
1853: This function is similar to `VecGetArrayRead()` which maps the local
1854: portion into a raw pointer. `VecGetLocalVectorRead()` is usually
1855: almost as efficient as `VecGetArrayRead()` but in certain circumstances
1856: `VecGetLocalVectorRead()` can be much more efficient than
1857: `VecGetArrayRead()`. This is because the construction of a contiguous
1858: array representing the vector data required by `VecGetArrayRead()` can
1859: be an expensive operation for certain vector types. For example, for
1860: GPU vectors `VecGetArrayRead()` requires that the data between device
1861: and host is synchronized.
1863: Unlike `VecGetLocalVector()`, this routine is not collective and
1864: preserves cached information.
1866: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1867: @*/
1868: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1869: {
1870: PetscFunctionBegin;
1873: VecCheckSameLocalSize(v, 1, w, 2);
1874: if (v->ops->getlocalvectorread) {
1875: PetscUseTypeMethod(v, getlocalvectorread, w);
1876: } else {
1877: PetscScalar *a;
1879: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1880: PetscCall(VecPlaceArray(w, a));
1881: }
1882: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1883: PetscCall(VecLockReadPush(v));
1884: PetscCall(VecLockReadPush(w));
1885: PetscFunctionReturn(PETSC_SUCCESS);
1886: }
1888: /*@
1889: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1890: previously mapped into a vector using `VecGetLocalVectorRead()`.
1892: Not Collective.
1894: Input Parameters:
1895: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
1896: - w - The vector into which the local portion of `v` was mapped.
1898: Level: beginner
1900: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1901: @*/
1902: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1903: {
1904: PetscFunctionBegin;
1907: if (v->ops->restorelocalvectorread) {
1908: PetscUseTypeMethod(v, restorelocalvectorread, w);
1909: } else {
1910: const PetscScalar *a;
1912: PetscCall(VecGetArrayRead(w, &a));
1913: PetscCall(VecRestoreArrayRead(v, &a));
1914: PetscCall(VecResetArray(w));
1915: }
1916: PetscCall(VecLockReadPop(v));
1917: PetscCall(VecLockReadPop(w));
1918: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1919: PetscFunctionReturn(PETSC_SUCCESS);
1920: }
1922: /*@
1923: VecGetLocalVector - Maps the local portion of a vector into a
1924: vector.
1926: Collective
1928: Input Parameter:
1929: . v - The vector for which the local vector is desired.
1931: Output Parameter:
1932: . w - Upon exit this contains the local vector.
1934: Level: beginner
1936: Notes:
1937: You must call `VecRestoreLocalVector()` when the local
1938: vector is no longer needed.
1940: This function is similar to `VecGetArray()` which maps the local
1941: portion into a raw pointer. `VecGetLocalVector()` is usually about as
1942: efficient as `VecGetArray()` but in certain circumstances
1943: `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1944: This is because the construction of a contiguous array representing
1945: the vector data required by `VecGetArray()` can be an expensive
1946: operation for certain vector types. For example, for GPU vectors
1947: `VecGetArray()` requires that the data between device and host is
1948: synchronized.
1950: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1951: @*/
1952: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1953: {
1954: PetscFunctionBegin;
1957: VecCheckSameLocalSize(v, 1, w, 2);
1958: if (v->ops->getlocalvector) {
1959: PetscUseTypeMethod(v, getlocalvector, w);
1960: } else {
1961: PetscScalar *a;
1963: PetscCall(VecGetArray(v, &a));
1964: PetscCall(VecPlaceArray(w, a));
1965: }
1966: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1967: PetscFunctionReturn(PETSC_SUCCESS);
1968: }
1970: /*@
1971: VecRestoreLocalVector - Unmaps the local portion of a vector
1972: previously mapped into a vector using `VecGetLocalVector()`.
1974: Logically Collective.
1976: Input Parameters:
1977: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
1978: - w - The vector into which the local portion of `v` was mapped.
1980: Level: beginner
1982: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1983: @*/
1984: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1985: {
1986: PetscFunctionBegin;
1989: if (v->ops->restorelocalvector) {
1990: PetscUseTypeMethod(v, restorelocalvector, w);
1991: } else {
1992: PetscScalar *a;
1993: PetscCall(VecGetArray(w, &a));
1994: PetscCall(VecRestoreArray(v, &a));
1995: PetscCall(VecResetArray(w));
1996: }
1997: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1998: PetscCall(PetscObjectStateIncrease((PetscObject)v));
1999: PetscFunctionReturn(PETSC_SUCCESS);
2000: }
2002: /*@C
2003: VecGetArray - Returns a pointer to a contiguous array that contains this
2004: MPI processes's portion of the vector data
2006: Logically Collective
2008: Input Parameter:
2009: . x - the vector
2011: Output Parameter:
2012: . a - location to put pointer to the array
2014: Level: beginner
2016: Notes:
2017: For the standard PETSc vectors, `VecGetArray()` returns a pointer to the local data array and
2018: does not use any copies. If the underlying vector data is not stored in a contiguous array
2019: this routine will copy the data to a contiguous array and return a pointer to that. You MUST
2020: call `VecRestoreArray()` when you no longer need access to the array.
2022: For vectors that may also have the array data in GPU memory, for example, `VECCUDA`, this call ensures the CPU array has the
2023: most recent array values by copying the data from the GPU memory if needed.
2025: Fortran Note:
2026: .vb
2027: PetscScalar, pointer :: a(:)
2028: .ve
2030: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecPlaceArray()`, `VecGetArray2d()`,
2031: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`, `VecGetArrayAndMemType()`
2032: @*/
2033: PetscErrorCode VecGetArray(Vec x, PetscScalar *a[])
2034: {
2035: PetscFunctionBegin;
2037: PetscCall(VecSetErrorIfLocked(x, 1));
2038: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
2039: PetscUseTypeMethod(x, getarray, a);
2040: } else if (x->petscnative) { /* VECSTANDARD */
2041: *a = *((PetscScalar **)x->data);
2042: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
2043: PetscFunctionReturn(PETSC_SUCCESS);
2044: }
2046: /*@C
2047: VecRestoreArray - Restores a vector after `VecGetArray()` has been called and the array is no longer needed
2049: Logically Collective
2051: Input Parameters:
2052: + x - the vector
2053: - a - location of pointer to array obtained from `VecGetArray()`
2055: Level: beginner
2057: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2058: `VecGetArrayPair()`, `VecRestoreArrayPair()`
2059: @*/
2060: PetscErrorCode VecRestoreArray(Vec x, PetscScalar *a[])
2061: {
2062: PetscFunctionBegin;
2064: if (a) PetscAssertPointer(a, 2);
2065: if (x->ops->restorearray) {
2066: PetscUseTypeMethod(x, restorearray, a);
2067: } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
2068: if (a) *a = NULL;
2069: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2070: PetscFunctionReturn(PETSC_SUCCESS);
2071: }
2072: /*@C
2073: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
2075: Not Collective
2077: Input Parameter:
2078: . x - the vector
2080: Output Parameter:
2081: . a - the array
2083: Level: beginner
2085: Notes:
2086: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2088: Unlike `VecGetArray()`, preserves cached information like vector norms.
2090: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
2091: implementations may require a copy, but such implementations should cache the contiguous representation so that
2092: only one copy is performed when this routine is called multiple times in sequence.
2094: For vectors that may also have the array data in GPU memory, for example, `VECCUDA`, this call ensures the CPU array has the
2095: most recent array values by copying the data from the GPU memory if needed.
2097: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2098: `VecGetArrayAndMemType()`
2099: @*/
2100: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar *a[])
2101: {
2102: PetscFunctionBegin;
2104: PetscAssertPointer(a, 2);
2105: if (x->ops->getarrayread) {
2106: PetscUseTypeMethod(x, getarrayread, a);
2107: } else if (x->ops->getarray) {
2108: PetscObjectState state;
2110: /* VECNEST, VECCUDA, VECKOKKOS etc */
2111: // x->ops->getarray may bump the object state, but since we know this is a read-only get
2112: // we can just undo that
2113: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2114: PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
2115: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2116: } else if (x->petscnative) {
2117: /* VECSTANDARD */
2118: *a = *((PetscScalar **)x->data);
2119: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
2120: PetscFunctionReturn(PETSC_SUCCESS);
2121: }
2123: /*@C
2124: VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`
2126: Not Collective
2128: Input Parameters:
2129: + x - the vector
2130: - a - the array
2132: Level: beginner
2134: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2135: @*/
2136: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar *a[])
2137: {
2138: PetscFunctionBegin;
2140: if (a) PetscAssertPointer(a, 2);
2141: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
2142: /* nothing */
2143: } else if (x->ops->restorearrayread) { /* VECNEST */
2144: PetscUseTypeMethod(x, restorearrayread, a);
2145: } else { /* No one? */
2146: PetscObjectState state;
2148: // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2149: // we can just undo that
2150: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2151: PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2152: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2153: }
2154: if (a) *a = NULL;
2155: PetscFunctionReturn(PETSC_SUCCESS);
2156: }
2158: /*@C
2159: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
2160: MPI processes's portion of the vector data.
2162: Logically Collective
2164: Input Parameter:
2165: . x - the vector
2167: Output Parameter:
2168: . a - location to put pointer to the array
2170: Level: intermediate
2172: Note:
2173: The values in this array are NOT valid, the caller of this routine is responsible for putting
2174: values into the array; any values it does not set will be invalid.
2176: The array must be returned using a matching call to `VecRestoreArrayWrite()`.
2178: For vectors associated with GPUs, the host and device vectors are not synchronized before
2179: giving access. If you need correct values in the array use `VecGetArray()`
2181: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecPlaceArray()`, `VecGetArray2d()`,
2182: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`, `VecGetArrayAndMemType()`
2183: @*/
2184: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar *a[])
2185: {
2186: PetscFunctionBegin;
2188: PetscAssertPointer(a, 2);
2189: PetscCall(VecSetErrorIfLocked(x, 1));
2190: if (x->ops->getarraywrite) {
2191: PetscUseTypeMethod(x, getarraywrite, a);
2192: } else {
2193: PetscCall(VecGetArray(x, a));
2194: }
2195: PetscFunctionReturn(PETSC_SUCCESS);
2196: }
2198: /*@C
2199: VecRestoreArrayWrite - Restores a vector after `VecGetArrayWrite()` has been called.
2201: Logically Collective
2203: Input Parameters:
2204: + x - the vector
2205: - a - location of pointer to array obtained from `VecGetArray()`
2207: Level: beginner
2209: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2210: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2211: @*/
2212: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar *a[])
2213: {
2214: PetscFunctionBegin;
2216: if (a) PetscAssertPointer(a, 2);
2217: if (x->ops->restorearraywrite) {
2218: PetscUseTypeMethod(x, restorearraywrite, a);
2219: } else if (x->ops->restorearray) {
2220: PetscUseTypeMethod(x, restorearray, a);
2221: }
2222: if (a) *a = NULL;
2223: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2224: PetscFunctionReturn(PETSC_SUCCESS);
2225: }
2227: /*@C
2228: VecGetArrays - Returns a pointer to the arrays in a set of vectors
2229: that were created by a call to `VecDuplicateVecs()`.
2231: Logically Collective; No Fortran Support
2233: Input Parameters:
2234: + x - the vectors
2235: - n - the number of vectors
2237: Output Parameter:
2238: . a - location to put pointer to the array
2240: Level: intermediate
2242: Note:
2243: You MUST call `VecRestoreArrays()` when you no longer need access to the arrays.
2245: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2246: @*/
2247: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2248: {
2249: PetscInt i;
2250: PetscScalar **q;
2252: PetscFunctionBegin;
2253: PetscAssertPointer(x, 1);
2255: PetscAssertPointer(a, 3);
2256: PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2257: PetscCall(PetscMalloc1(n, &q));
2258: for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2259: *a = q;
2260: PetscFunctionReturn(PETSC_SUCCESS);
2261: }
2263: /*@C
2264: VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2265: has been called.
2267: Logically Collective; No Fortran Support
2269: Input Parameters:
2270: + x - the vector
2271: . n - the number of vectors
2272: - a - location of pointer to arrays obtained from `VecGetArrays()`
2274: Notes:
2275: For regular PETSc vectors this routine does not involve any copies. For
2276: any special vectors that do not store local vector data in a contiguous
2277: array, this routine will copy the data back into the underlying
2278: vector data structure from the arrays obtained with `VecGetArrays()`.
2280: Level: intermediate
2282: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2283: @*/
2284: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2285: {
2286: PetscInt i;
2287: PetscScalar **q = *a;
2289: PetscFunctionBegin;
2290: PetscAssertPointer(x, 1);
2292: PetscAssertPointer(a, 3);
2294: for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2295: PetscCall(PetscFree(q));
2296: PetscFunctionReturn(PETSC_SUCCESS);
2297: }
2299: /*@C
2300: VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g.,
2301: `VECCUDA`), the returned pointer will be a device pointer to the device memory that contains
2302: this MPI processes's portion of the vector data.
2304: Logically Collective; No Fortran Support
2306: Input Parameter:
2307: . x - the vector
2309: Output Parameters:
2310: + a - location to put pointer to the array
2311: - mtype - memory type of the array
2313: Level: beginner
2315: Note:
2316: Device data is guaranteed to have the latest value. Otherwise, when this is a host vector
2317: (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host
2318: pointer.
2320: For `VECKOKKOS`, if Kokkos is configured without device (e.g., use serial or openmp), per
2321: this function, the vector works like `VECSEQ`/`VECMPI`; otherwise, it works like `VECCUDA` or
2322: `VECHIP` etc.
2324: Use `VecRestoreArrayAndMemType()` when the array access is no longer needed.
2326: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`,
2327: `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2328: @*/
2329: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar *a[], PetscMemType *mtype)
2330: {
2331: PetscFunctionBegin;
2334: if (a) PetscAssertPointer(a, 2);
2335: if (mtype) PetscAssertPointer(mtype, 3);
2336: PetscCall(VecSetErrorIfLocked(x, 1));
2337: if (x->ops->getarrayandmemtype) {
2338: /* VECCUDA, VECKOKKOS etc */
2339: PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2340: } else {
2341: /* VECSTANDARD, VECNEST, VECVIENNACL */
2342: PetscCall(VecGetArray(x, a));
2343: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2344: }
2345: PetscFunctionReturn(PETSC_SUCCESS);
2346: }
2348: /*@C
2349: VecRestoreArrayAndMemType - Restores a vector after `VecGetArrayAndMemType()` has been called.
2351: Logically Collective; No Fortran Support
2353: Input Parameters:
2354: + x - the vector
2355: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`
2357: Level: beginner
2359: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`,
2360: `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2361: @*/
2362: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar *a[])
2363: {
2364: PetscFunctionBegin;
2367: if (a) PetscAssertPointer(a, 2);
2368: if (x->ops->restorearrayandmemtype) {
2369: /* VECCUDA, VECKOKKOS etc */
2370: PetscUseTypeMethod(x, restorearrayandmemtype, a);
2371: } else {
2372: /* VECNEST, VECVIENNACL */
2373: PetscCall(VecRestoreArray(x, a));
2374: } /* VECSTANDARD does nothing */
2375: if (a) *a = NULL;
2376: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2377: PetscFunctionReturn(PETSC_SUCCESS);
2378: }
2380: /*@C
2381: VecGetArrayReadAndMemType - Like `VecGetArrayRead()`, but if the input vector is a device vector, it will return a read-only device pointer.
2382: The returned pointer is guaranteed to point to up-to-date data. For host vectors, it functions as `VecGetArrayRead()`.
2384: Not Collective; No Fortran Support
2386: Input Parameter:
2387: . x - the vector
2389: Output Parameters:
2390: + a - the array
2391: - mtype - memory type of the array
2393: Level: beginner
2395: Notes:
2396: The array must be returned using a matching call to `VecRestoreArrayReadAndMemType()`.
2398: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2399: @*/
2400: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar *a[], PetscMemType *mtype)
2401: {
2402: PetscFunctionBegin;
2405: PetscAssertPointer(a, 2);
2406: if (mtype) PetscAssertPointer(mtype, 3);
2407: if (x->ops->getarrayreadandmemtype) {
2408: /* VECCUDA/VECHIP though they are also petscnative */
2409: PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2410: } else if (x->ops->getarrayandmemtype) {
2411: /* VECKOKKOS */
2412: PetscObjectState state;
2414: // see VecGetArrayRead() for why
2415: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2416: PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2417: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2418: } else {
2419: PetscCall(VecGetArrayRead(x, a));
2420: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2421: }
2422: PetscFunctionReturn(PETSC_SUCCESS);
2423: }
2425: /*@C
2426: VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`
2428: Not Collective; No Fortran Support
2430: Input Parameters:
2431: + x - the vector
2432: - a - the array
2434: Level: beginner
2436: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2437: @*/
2438: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar *a[])
2439: {
2440: PetscFunctionBegin;
2443: if (a) PetscAssertPointer(a, 2);
2444: if (x->ops->restorearrayreadandmemtype) {
2445: /* VECCUDA/VECHIP */
2446: PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2447: } else if (!x->petscnative) {
2448: /* VECNEST */
2449: PetscCall(VecRestoreArrayRead(x, a));
2450: }
2451: if (a) *a = NULL;
2452: PetscFunctionReturn(PETSC_SUCCESS);
2453: }
2455: /*@C
2456: VecGetArrayWriteAndMemType - Like `VecGetArrayWrite()`, but if this is a device vector it will always return
2457: a device pointer to the device memory that contains this processor's portion of the vector data.
2459: Logically Collective; No Fortran Support
2461: Input Parameter:
2462: . x - the vector
2464: Output Parameters:
2465: + a - the array
2466: - mtype - memory type of the array
2468: Level: beginner
2470: Note:
2471: The array must be returned using a matching call to `VecRestoreArrayWriteAndMemType()`, where it will label the device memory as most recent.
2473: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2474: @*/
2475: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar *a[], PetscMemType *mtype)
2476: {
2477: PetscFunctionBegin;
2480: PetscCall(VecSetErrorIfLocked(x, 1));
2481: PetscAssertPointer(a, 2);
2482: if (mtype) PetscAssertPointer(mtype, 3);
2483: if (x->ops->getarraywriteandmemtype) {
2484: /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2485: PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2486: } else if (x->ops->getarrayandmemtype) {
2487: PetscCall(VecGetArrayAndMemType(x, a, mtype));
2488: } else {
2489: /* VECNEST, VECVIENNACL */
2490: PetscCall(VecGetArrayWrite(x, a));
2491: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2492: }
2493: PetscFunctionReturn(PETSC_SUCCESS);
2494: }
2496: /*@C
2497: VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`
2499: Logically Collective; No Fortran Support
2501: Input Parameters:
2502: + x - the vector
2503: - a - the array
2505: Level: beginner
2507: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2508: @*/
2509: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar *a[])
2510: {
2511: PetscFunctionBegin;
2514: PetscCall(VecSetErrorIfLocked(x, 1));
2515: if (a) PetscAssertPointer(a, 2);
2516: if (x->ops->restorearraywriteandmemtype) {
2517: /* VECCUDA/VECHIP */
2518: PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2519: PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2520: } else if (x->ops->restorearrayandmemtype) {
2521: PetscCall(VecRestoreArrayAndMemType(x, a));
2522: } else {
2523: PetscCall(VecRestoreArray(x, a));
2524: }
2525: if (a) *a = NULL;
2526: PetscFunctionReturn(PETSC_SUCCESS);
2527: }
2529: /*@
2530: VecPlaceArray - Allows one to replace the array in a vector with an
2531: array provided by the user. This is useful to avoid copying an array
2532: into a vector.
2534: Logically Collective
2536: Input Parameters:
2537: + vec - the vector
2538: - array - the array
2540: Level: developer
2542: Notes:
2543: Adding `const` to `array` was an oversight, as subsequent operations on `vec` would
2544: likely modify the data in `array`. However, we have kept it to avoid breaking APIs.
2546: Use `VecReplaceArray()` instead to permanently replace the array
2548: You can return to the original array with a call to `VecResetArray()`. `vec` does not take
2549: ownership of `array` in any way.
2551: The user must free `array` themselves but be careful not to
2552: do so before the vector has either been destroyed, had its original array restored with
2553: `VecResetArray()` or permanently replaced with `VecReplaceArray()`.
2555: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2556: @*/
2557: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2558: {
2559: PetscFunctionBegin;
2562: if (array) PetscAssertPointer(array, 2);
2563: PetscUseTypeMethod(vec, placearray, array);
2564: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2565: PetscFunctionReturn(PETSC_SUCCESS);
2566: }
2568: /*@C
2569: VecReplaceArray - Allows one to replace the array in a vector with an
2570: array provided by the user. This is useful to avoid copying an array
2571: into a vector.
2573: Logically Collective; No Fortran Support
2575: Input Parameters:
2576: + vec - the vector
2577: - array - the array
2579: Level: developer
2581: Notes:
2582: Adding `const` to `array` was an oversight, as subsequent operations on `vec` would
2583: likely modify the data in `array`. However, we have kept it to avoid breaking APIs.
2585: This permanently replaces the array and frees the memory associated
2586: with the old array. Use `VecPlaceArray()` to temporarily replace the array.
2588: The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
2589: freed by the user. It will be freed when the vector is destroyed.
2591: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2592: @*/
2593: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2594: {
2595: PetscFunctionBegin;
2598: PetscUseTypeMethod(vec, replacearray, array);
2599: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2600: PetscFunctionReturn(PETSC_SUCCESS);
2601: }
2603: /*@C
2604: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2605: processor's portion of the vector data. You MUST call `VecRestoreArray2d()`
2606: when you no longer need access to the array.
2608: Logically Collective
2610: Input Parameters:
2611: + x - the vector
2612: . m - first dimension of two dimensional array
2613: . n - second dimension of two dimensional array
2614: . mstart - first index you will use in first coordinate direction (often 0)
2615: - nstart - first index in the second coordinate direction (often 0)
2617: Output Parameter:
2618: . a - location to put pointer to the array
2620: Level: developer
2622: Notes:
2623: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2624: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2625: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2626: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2628: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2630: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2631: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2632: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2633: @*/
2634: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2635: {
2636: PetscInt i, N;
2637: PetscScalar *aa;
2639: PetscFunctionBegin;
2641: PetscAssertPointer(a, 6);
2643: PetscCall(VecGetLocalSize(x, &N));
2644: 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);
2645: PetscCall(VecGetArray(x, &aa));
2647: PetscCall(PetscMalloc1(m, a));
2648: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2649: *a -= mstart;
2650: PetscFunctionReturn(PETSC_SUCCESS);
2651: }
2653: /*@C
2654: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2655: processor's portion of the vector data. You MUST call `VecRestoreArray2dWrite()`
2656: when you no longer need access to the array.
2658: Logically Collective
2660: Input Parameters:
2661: + x - the vector
2662: . m - first dimension of two dimensional array
2663: . n - second dimension of two dimensional array
2664: . mstart - first index you will use in first coordinate direction (often 0)
2665: - nstart - first index in the second coordinate direction (often 0)
2667: Output Parameter:
2668: . a - location to put pointer to the array
2670: Level: developer
2672: Notes:
2673: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2674: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2675: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2676: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2678: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2680: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2681: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2682: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2683: @*/
2684: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2685: {
2686: PetscInt i, N;
2687: PetscScalar *aa;
2689: PetscFunctionBegin;
2691: PetscAssertPointer(a, 6);
2693: PetscCall(VecGetLocalSize(x, &N));
2694: 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);
2695: PetscCall(VecGetArrayWrite(x, &aa));
2697: PetscCall(PetscMalloc1(m, a));
2698: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2699: *a -= mstart;
2700: PetscFunctionReturn(PETSC_SUCCESS);
2701: }
2703: /*@C
2704: VecRestoreArray2d - Restores a vector after `VecGetArray2d()` has been called.
2706: Logically Collective
2708: Input Parameters:
2709: + x - the vector
2710: . m - first dimension of two dimensional array
2711: . n - second dimension of the two dimensional array
2712: . mstart - first index you will use in first coordinate direction (often 0)
2713: . nstart - first index in the second coordinate direction (often 0)
2714: - a - location of pointer to array obtained from `VecGetArray2d()`
2716: Level: developer
2718: Notes:
2719: For regular PETSc vectors this routine does not involve any copies. For
2720: any special vectors that do not store local vector data in a contiguous
2721: array, this routine will copy the data back into the underlying
2722: vector data structure from the array obtained with `VecGetArray()`.
2724: This routine actually zeros out the `a` pointer.
2726: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2727: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
2728: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2729: @*/
2730: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2731: {
2732: void *dummy;
2734: PetscFunctionBegin;
2736: PetscAssertPointer(a, 6);
2738: dummy = (void *)(*a + mstart);
2739: PetscCall(PetscFree(dummy));
2740: PetscCall(VecRestoreArray(x, NULL));
2741: *a = NULL;
2742: PetscFunctionReturn(PETSC_SUCCESS);
2743: }
2745: /*@C
2746: VecRestoreArray2dWrite - Restores a vector after `VecGetArray2dWrite()` has been called.
2748: Logically Collective
2750: Input Parameters:
2751: + x - the vector
2752: . m - first dimension of two dimensional array
2753: . n - second dimension of the two dimensional array
2754: . mstart - first index you will use in first coordinate direction (often 0)
2755: . nstart - first index in the second coordinate direction (often 0)
2756: - a - location of pointer to array obtained from `VecGetArray2d()`
2758: Level: developer
2760: Notes:
2761: For regular PETSc vectors this routine does not involve any copies. For
2762: any special vectors that do not store local vector data in a contiguous
2763: array, this routine will copy the data back into the underlying
2764: vector data structure from the array obtained with `VecGetArray()`.
2766: This routine actually zeros out the `a` pointer.
2768: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2769: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
2770: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2771: @*/
2772: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2773: {
2774: void *dummy;
2776: PetscFunctionBegin;
2778: PetscAssertPointer(a, 6);
2780: dummy = (void *)(*a + mstart);
2781: PetscCall(PetscFree(dummy));
2782: PetscCall(VecRestoreArrayWrite(x, NULL));
2783: PetscFunctionReturn(PETSC_SUCCESS);
2784: }
2786: /*@C
2787: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2788: processor's portion of the vector data. You MUST call `VecRestoreArray1d()`
2789: when you no longer need access to the array.
2791: Logically Collective
2793: Input Parameters:
2794: + x - the vector
2795: . m - first dimension of two dimensional array
2796: - mstart - first index you will use in first coordinate direction (often 0)
2798: Output Parameter:
2799: . a - location to put pointer to the array
2801: Level: developer
2803: Notes:
2804: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2805: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2806: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2808: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2810: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2811: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2812: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2813: @*/
2814: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2815: {
2816: PetscInt N;
2818: PetscFunctionBegin;
2820: PetscAssertPointer(a, 4);
2822: PetscCall(VecGetLocalSize(x, &N));
2823: 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);
2824: PetscCall(VecGetArray(x, a));
2825: *a -= mstart;
2826: PetscFunctionReturn(PETSC_SUCCESS);
2827: }
2829: /*@C
2830: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
2831: processor's portion of the vector data. You MUST call `VecRestoreArray1dWrite()`
2832: when you no longer need access to the array.
2834: Logically Collective
2836: Input Parameters:
2837: + x - the vector
2838: . m - first dimension of two dimensional array
2839: - mstart - first index you will use in first coordinate direction (often 0)
2841: Output Parameter:
2842: . a - location to put pointer to the array
2844: Level: developer
2846: Notes:
2847: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2848: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2849: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2851: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2853: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2854: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2855: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2856: @*/
2857: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2858: {
2859: PetscInt N;
2861: PetscFunctionBegin;
2863: PetscAssertPointer(a, 4);
2865: PetscCall(VecGetLocalSize(x, &N));
2866: 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);
2867: PetscCall(VecGetArrayWrite(x, a));
2868: *a -= mstart;
2869: PetscFunctionReturn(PETSC_SUCCESS);
2870: }
2872: /*@C
2873: VecRestoreArray1d - Restores a vector after `VecGetArray1d()` has been called.
2875: Logically Collective
2877: Input Parameters:
2878: + x - the vector
2879: . m - first dimension of two dimensional array
2880: . mstart - first index you will use in first coordinate direction (often 0)
2881: - a - location of pointer to array obtained from `VecGetArray1d()`
2883: Level: developer
2885: Notes:
2886: For regular PETSc vectors this routine does not involve any copies. For
2887: any special vectors that do not store local vector data in a contiguous
2888: array, this routine will copy the data back into the underlying
2889: vector data structure from the array obtained with `VecGetArray1d()`.
2891: This routine actually zeros out the `a` pointer.
2893: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2894: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
2895: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2896: @*/
2897: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2898: {
2899: PetscFunctionBegin;
2902: PetscCall(VecRestoreArray(x, NULL));
2903: *a = NULL;
2904: PetscFunctionReturn(PETSC_SUCCESS);
2905: }
2907: /*@C
2908: VecRestoreArray1dWrite - Restores a vector after `VecGetArray1dWrite()` has been called.
2910: Logically Collective
2912: Input Parameters:
2913: + x - the vector
2914: . m - first dimension of two dimensional array
2915: . mstart - first index you will use in first coordinate direction (often 0)
2916: - a - location of pointer to array obtained from `VecGetArray1d()`
2918: Level: developer
2920: Notes:
2921: For regular PETSc vectors this routine does not involve any copies. For
2922: any special vectors that do not store local vector data in a contiguous
2923: array, this routine will copy the data back into the underlying
2924: vector data structure from the array obtained with `VecGetArray1d()`.
2926: This routine actually zeros out the `a` pointer.
2928: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
2929: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
2930: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2931: @*/
2932: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2933: {
2934: PetscFunctionBegin;
2937: PetscCall(VecRestoreArrayWrite(x, NULL));
2938: *a = NULL;
2939: PetscFunctionReturn(PETSC_SUCCESS);
2940: }
2942: /*@C
2943: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
2944: processor's portion of the vector data. You MUST call `VecRestoreArray3d()`
2945: when you no longer need access to the array.
2947: Logically Collective
2949: Input Parameters:
2950: + x - the vector
2951: . m - first dimension of three dimensional array
2952: . n - second dimension of three dimensional array
2953: . p - third dimension of three dimensional array
2954: . mstart - first index you will use in first coordinate direction (often 0)
2955: . nstart - first index in the second coordinate direction (often 0)
2956: - pstart - first index in the third coordinate direction (often 0)
2958: Output Parameter:
2959: . a - location to put pointer to the array
2961: Level: developer
2963: Notes:
2964: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
2965: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2966: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2967: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
2969: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2971: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
2972: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
2973: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2974: @*/
2975: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
2976: {
2977: PetscInt i, N, j;
2978: PetscScalar *aa, **b;
2980: PetscFunctionBegin;
2982: PetscAssertPointer(a, 8);
2984: PetscCall(VecGetLocalSize(x, &N));
2985: 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);
2986: PetscCall(VecGetArray(x, &aa));
2988: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
2989: b = (PetscScalar **)((*a) + m);
2990: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
2991: for (i = 0; i < m; i++)
2992: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset(aa, i * n * p + j * p - pstart);
2993: *a -= mstart;
2994: PetscFunctionReturn(PETSC_SUCCESS);
2995: }
2997: /*@C
2998: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
2999: processor's portion of the vector data. You MUST call `VecRestoreArray3dWrite()`
3000: when you no longer need access to the array.
3002: Logically Collective
3004: Input Parameters:
3005: + x - the vector
3006: . m - first dimension of three dimensional array
3007: . n - second dimension of three dimensional array
3008: . p - third dimension of three dimensional array
3009: . mstart - first index you will use in first coordinate direction (often 0)
3010: . nstart - first index in the second coordinate direction (often 0)
3011: - pstart - first index in the third coordinate direction (often 0)
3013: Output Parameter:
3014: . a - location to put pointer to the array
3016: Level: developer
3018: Notes:
3019: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3020: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3021: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3022: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3024: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3026: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3027: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3028: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3029: @*/
3030: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3031: {
3032: PetscInt i, N, j;
3033: PetscScalar *aa, **b;
3035: PetscFunctionBegin;
3037: PetscAssertPointer(a, 8);
3039: PetscCall(VecGetLocalSize(x, &N));
3040: 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);
3041: PetscCall(VecGetArrayWrite(x, &aa));
3043: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3044: b = (PetscScalar **)((*a) + m);
3045: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3046: for (i = 0; i < m; i++)
3047: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3049: *a -= mstart;
3050: PetscFunctionReturn(PETSC_SUCCESS);
3051: }
3053: /*@C
3054: VecRestoreArray3d - Restores a vector after `VecGetArray3d()` has been called.
3056: Logically Collective
3058: Input Parameters:
3059: + x - the vector
3060: . m - first dimension of three dimensional array
3061: . n - second dimension of the three dimensional array
3062: . p - third dimension of the three dimensional array
3063: . mstart - first index you will use in first coordinate direction (often 0)
3064: . nstart - first index in the second coordinate direction (often 0)
3065: . pstart - first index in the third coordinate direction (often 0)
3066: - a - location of pointer to array obtained from VecGetArray3d()
3068: Level: developer
3070: Notes:
3071: For regular PETSc vectors this routine does not involve any copies. For
3072: any special vectors that do not store local vector data in a contiguous
3073: array, this routine will copy the data back into the underlying
3074: vector data structure from the array obtained with `VecGetArray()`.
3076: This routine actually zeros out the `a` pointer.
3078: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3079: `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3080: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3081: @*/
3082: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3083: {
3084: void *dummy;
3086: PetscFunctionBegin;
3088: PetscAssertPointer(a, 8);
3090: dummy = (void *)(*a + mstart);
3091: PetscCall(PetscFree(dummy));
3092: PetscCall(VecRestoreArray(x, NULL));
3093: *a = NULL;
3094: PetscFunctionReturn(PETSC_SUCCESS);
3095: }
3097: /*@C
3098: VecRestoreArray3dWrite - Restores a vector after `VecGetArray3dWrite()` has been called.
3100: Logically Collective
3102: Input Parameters:
3103: + x - the vector
3104: . m - first dimension of three dimensional array
3105: . n - second dimension of the three dimensional array
3106: . p - third dimension of the three dimensional array
3107: . mstart - first index you will use in first coordinate direction (often 0)
3108: . nstart - first index in the second coordinate direction (often 0)
3109: . pstart - first index in the third coordinate direction (often 0)
3110: - a - location of pointer to array obtained from VecGetArray3d()
3112: Level: developer
3114: Notes:
3115: For regular PETSc vectors this routine does not involve any copies. For
3116: any special vectors that do not store local vector data in a contiguous
3117: array, this routine will copy the data back into the underlying
3118: vector data structure from the array obtained with `VecGetArray()`.
3120: This routine actually zeros out the `a` pointer.
3122: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3123: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3124: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3125: @*/
3126: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3127: {
3128: void *dummy;
3130: PetscFunctionBegin;
3132: PetscAssertPointer(a, 8);
3134: dummy = (void *)(*a + mstart);
3135: PetscCall(PetscFree(dummy));
3136: PetscCall(VecRestoreArrayWrite(x, NULL));
3137: *a = NULL;
3138: PetscFunctionReturn(PETSC_SUCCESS);
3139: }
3141: /*@C
3142: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3143: processor's portion of the vector data. You MUST call `VecRestoreArray4d()`
3144: when you no longer need access to the array.
3146: Logically Collective
3148: Input Parameters:
3149: + x - the vector
3150: . m - first dimension of four dimensional array
3151: . n - second dimension of four dimensional array
3152: . p - third dimension of four dimensional array
3153: . q - fourth dimension of four dimensional array
3154: . mstart - first index you will use in first coordinate direction (often 0)
3155: . nstart - first index in the second coordinate direction (often 0)
3156: . pstart - first index in the third coordinate direction (often 0)
3157: - qstart - first index in the fourth coordinate direction (often 0)
3159: Output Parameter:
3160: . a - location to put pointer to the array
3162: Level: developer
3164: Notes:
3165: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3166: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3167: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3168: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3170: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3172: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3173: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3174: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3175: @*/
3176: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3177: {
3178: PetscInt i, N, j, k;
3179: PetscScalar *aa, ***b, **c;
3181: PetscFunctionBegin;
3183: PetscAssertPointer(a, 10);
3185: PetscCall(VecGetLocalSize(x, &N));
3186: 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);
3187: PetscCall(VecGetArray(x, &aa));
3189: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3190: b = (PetscScalar ***)((*a) + m);
3191: c = (PetscScalar **)(b + m * n);
3192: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3193: for (i = 0; i < m; i++)
3194: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3195: for (i = 0; i < m; i++)
3196: for (j = 0; j < n; j++)
3197: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3198: *a -= mstart;
3199: PetscFunctionReturn(PETSC_SUCCESS);
3200: }
3202: /*@C
3203: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3204: processor's portion of the vector data. You MUST call `VecRestoreArray4dWrite()`
3205: when you no longer need access to the array.
3207: Logically Collective
3209: Input Parameters:
3210: + x - the vector
3211: . m - first dimension of four dimensional array
3212: . n - second dimension of four dimensional array
3213: . p - third dimension of four dimensional array
3214: . q - fourth dimension of four dimensional array
3215: . mstart - first index you will use in first coordinate direction (often 0)
3216: . nstart - first index in the second coordinate direction (often 0)
3217: . pstart - first index in the third coordinate direction (often 0)
3218: - qstart - first index in the fourth coordinate direction (often 0)
3220: Output Parameter:
3221: . a - location to put pointer to the array
3223: Level: developer
3225: Notes:
3226: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3227: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3228: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3229: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3231: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3233: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3234: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3235: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3236: @*/
3237: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3238: {
3239: PetscInt i, N, j, k;
3240: PetscScalar *aa, ***b, **c;
3242: PetscFunctionBegin;
3244: PetscAssertPointer(a, 10);
3246: PetscCall(VecGetLocalSize(x, &N));
3247: 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);
3248: PetscCall(VecGetArrayWrite(x, &aa));
3250: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3251: b = (PetscScalar ***)((*a) + m);
3252: c = (PetscScalar **)(b + m * n);
3253: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3254: for (i = 0; i < m; i++)
3255: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3256: for (i = 0; i < m; i++)
3257: for (j = 0; j < n; j++)
3258: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3259: *a -= mstart;
3260: PetscFunctionReturn(PETSC_SUCCESS);
3261: }
3263: /*@C
3264: VecRestoreArray4d - Restores a vector after `VecGetArray4d()` has been called.
3266: Logically Collective
3268: Input Parameters:
3269: + x - the vector
3270: . m - first dimension of four dimensional array
3271: . n - second dimension of the four dimensional array
3272: . p - third dimension of the four dimensional array
3273: . q - fourth dimension of the four dimensional array
3274: . mstart - first index you will use in first coordinate direction (often 0)
3275: . nstart - first index in the second coordinate direction (often 0)
3276: . pstart - first index in the third coordinate direction (often 0)
3277: . qstart - first index in the fourth coordinate direction (often 0)
3278: - a - location of pointer to array obtained from VecGetArray4d()
3280: Level: developer
3282: Notes:
3283: For regular PETSc vectors this routine does not involve any copies. For
3284: any special vectors that do not store local vector data in a contiguous
3285: array, this routine will copy the data back into the underlying
3286: vector data structure from the array obtained with `VecGetArray()`.
3288: This routine actually zeros out the `a` pointer.
3290: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3291: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3292: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`
3293: @*/
3294: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3295: {
3296: void *dummy;
3298: PetscFunctionBegin;
3300: PetscAssertPointer(a, 10);
3302: dummy = (void *)(*a + mstart);
3303: PetscCall(PetscFree(dummy));
3304: PetscCall(VecRestoreArray(x, NULL));
3305: *a = NULL;
3306: PetscFunctionReturn(PETSC_SUCCESS);
3307: }
3309: /*@C
3310: VecRestoreArray4dWrite - Restores a vector after `VecGetArray4dWrite()` has been called.
3312: Logically Collective
3314: Input Parameters:
3315: + x - the vector
3316: . m - first dimension of four dimensional array
3317: . n - second dimension of the four dimensional array
3318: . p - third dimension of the four dimensional array
3319: . q - fourth dimension of the four dimensional array
3320: . mstart - first index you will use in first coordinate direction (often 0)
3321: . nstart - first index in the second coordinate direction (often 0)
3322: . pstart - first index in the third coordinate direction (often 0)
3323: . qstart - first index in the fourth coordinate direction (often 0)
3324: - a - location of pointer to array obtained from `VecGetArray4d()`
3326: Level: developer
3328: Notes:
3329: For regular PETSc vectors this routine does not involve any copies. For
3330: any special vectors that do not store local vector data in a contiguous
3331: array, this routine will copy the data back into the underlying
3332: vector data structure from the array obtained with `VecGetArray()`.
3334: This routine actually zeros out the `a` pointer.
3336: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3337: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3338: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3339: @*/
3340: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3341: {
3342: void *dummy;
3344: PetscFunctionBegin;
3346: PetscAssertPointer(a, 10);
3348: dummy = (void *)(*a + mstart);
3349: PetscCall(PetscFree(dummy));
3350: PetscCall(VecRestoreArrayWrite(x, NULL));
3351: *a = NULL;
3352: PetscFunctionReturn(PETSC_SUCCESS);
3353: }
3355: /*@C
3356: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3357: processor's portion of the vector data. You MUST call `VecRestoreArray2dRead()`
3358: when you no longer need access to the array.
3360: Logically Collective
3362: Input Parameters:
3363: + x - the vector
3364: . m - first dimension of two dimensional array
3365: . n - second dimension of two dimensional array
3366: . mstart - first index you will use in first coordinate direction (often 0)
3367: - nstart - first index in the second coordinate direction (often 0)
3369: Output Parameter:
3370: . a - location to put pointer to the array
3372: Level: developer
3374: Notes:
3375: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
3376: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3377: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3378: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
3380: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3382: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3383: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3384: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3385: @*/
3386: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3387: {
3388: PetscInt i, N;
3389: const PetscScalar *aa;
3391: PetscFunctionBegin;
3393: PetscAssertPointer(a, 6);
3395: PetscCall(VecGetLocalSize(x, &N));
3396: 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);
3397: PetscCall(VecGetArrayRead(x, &aa));
3399: PetscCall(PetscMalloc1(m, a));
3400: for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3401: *a -= mstart;
3402: PetscFunctionReturn(PETSC_SUCCESS);
3403: }
3405: /*@C
3406: VecRestoreArray2dRead - Restores a vector after `VecGetArray2dRead()` has been called.
3408: Logically Collective
3410: Input Parameters:
3411: + x - the vector
3412: . m - first dimension of two dimensional array
3413: . n - second dimension of the two dimensional array
3414: . mstart - first index you will use in first coordinate direction (often 0)
3415: . nstart - first index in the second coordinate direction (often 0)
3416: - a - location of pointer to array obtained from VecGetArray2d()
3418: Level: developer
3420: Notes:
3421: For regular PETSc vectors this routine does not involve any copies. For
3422: any special vectors that do not store local vector data in a contiguous
3423: array, this routine will copy the data back into the underlying
3424: vector data structure from the array obtained with `VecGetArray()`.
3426: This routine actually zeros out the `a` pointer.
3428: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3429: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3430: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3431: @*/
3432: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3433: {
3434: void *dummy;
3436: PetscFunctionBegin;
3438: PetscAssertPointer(a, 6);
3440: dummy = (void *)(*a + mstart);
3441: PetscCall(PetscFree(dummy));
3442: PetscCall(VecRestoreArrayRead(x, NULL));
3443: *a = NULL;
3444: PetscFunctionReturn(PETSC_SUCCESS);
3445: }
3447: /*@C
3448: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3449: processor's portion of the vector data. You MUST call `VecRestoreArray1dRead()`
3450: when you no longer need access to the array.
3452: Logically Collective
3454: Input Parameters:
3455: + x - the vector
3456: . m - first dimension of two dimensional array
3457: - mstart - first index you will use in first coordinate direction (often 0)
3459: Output Parameter:
3460: . a - location to put pointer to the array
3462: Level: developer
3464: Notes:
3465: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3466: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3467: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3469: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3471: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3472: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3473: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3474: @*/
3475: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3476: {
3477: PetscInt N;
3479: PetscFunctionBegin;
3481: PetscAssertPointer(a, 4);
3483: PetscCall(VecGetLocalSize(x, &N));
3484: 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);
3485: PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3486: *a -= mstart;
3487: PetscFunctionReturn(PETSC_SUCCESS);
3488: }
3490: /*@C
3491: VecRestoreArray1dRead - Restores a vector after `VecGetArray1dRead()` has been called.
3493: Logically Collective
3495: Input Parameters:
3496: + x - the vector
3497: . m - first dimension of two dimensional array
3498: . mstart - first index you will use in first coordinate direction (often 0)
3499: - a - location of pointer to array obtained from `VecGetArray1dRead()`
3501: Level: developer
3503: Notes:
3504: For regular PETSc vectors this routine does not involve any copies. For
3505: any special vectors that do not store local vector data in a contiguous
3506: array, this routine will copy the data back into the underlying
3507: vector data structure from the array obtained with `VecGetArray1dRead()`.
3509: This routine actually zeros out the `a` pointer.
3511: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3512: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3513: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3514: @*/
3515: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3516: {
3517: PetscFunctionBegin;
3520: PetscCall(VecRestoreArrayRead(x, NULL));
3521: *a = NULL;
3522: PetscFunctionReturn(PETSC_SUCCESS);
3523: }
3525: /*@C
3526: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3527: processor's portion of the vector data. You MUST call `VecRestoreArray3dRead()`
3528: when you no longer need access to the array.
3530: Logically Collective
3532: Input Parameters:
3533: + x - the vector
3534: . m - first dimension of three dimensional array
3535: . n - second dimension of three dimensional array
3536: . p - third dimension of three dimensional array
3537: . mstart - first index you will use in first coordinate direction (often 0)
3538: . nstart - first index in the second coordinate direction (often 0)
3539: - pstart - first index in the third coordinate direction (often 0)
3541: Output Parameter:
3542: . a - location to put pointer to the array
3544: Level: developer
3546: Notes:
3547: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3548: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3549: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3550: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3dRead()`.
3552: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3554: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3555: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3556: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3557: @*/
3558: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3559: {
3560: PetscInt i, N, j;
3561: const PetscScalar *aa;
3562: PetscScalar **b;
3564: PetscFunctionBegin;
3566: PetscAssertPointer(a, 8);
3568: PetscCall(VecGetLocalSize(x, &N));
3569: 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);
3570: PetscCall(VecGetArrayRead(x, &aa));
3572: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3573: b = (PetscScalar **)((*a) + m);
3574: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3575: for (i = 0; i < m; i++)
3576: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset((PetscScalar *)aa, i * n * p + j * p - pstart);
3577: *a -= mstart;
3578: PetscFunctionReturn(PETSC_SUCCESS);
3579: }
3581: /*@C
3582: VecRestoreArray3dRead - Restores a vector after `VecGetArray3dRead()` has been called.
3584: Logically Collective
3586: Input Parameters:
3587: + x - the vector
3588: . m - first dimension of three dimensional array
3589: . n - second dimension of the three dimensional array
3590: . p - third dimension of the three dimensional array
3591: . mstart - first index you will use in first coordinate direction (often 0)
3592: . nstart - first index in the second coordinate direction (often 0)
3593: . pstart - first index in the third coordinate direction (often 0)
3594: - a - location of pointer to array obtained from `VecGetArray3dRead()`
3596: Level: developer
3598: Notes:
3599: For regular PETSc vectors this routine does not involve any copies. For
3600: any special vectors that do not store local vector data in a contiguous
3601: array, this routine will copy the data back into the underlying
3602: vector data structure from the array obtained with `VecGetArray()`.
3604: This routine actually zeros out the `a` pointer.
3606: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3607: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3608: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3609: @*/
3610: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3611: {
3612: void *dummy;
3614: PetscFunctionBegin;
3616: PetscAssertPointer(a, 8);
3618: dummy = (void *)(*a + mstart);
3619: PetscCall(PetscFree(dummy));
3620: PetscCall(VecRestoreArrayRead(x, NULL));
3621: *a = NULL;
3622: PetscFunctionReturn(PETSC_SUCCESS);
3623: }
3625: /*@C
3626: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3627: processor's portion of the vector data. You MUST call `VecRestoreArray4dRead()`
3628: when you no longer need access to the array.
3630: Logically Collective
3632: Input Parameters:
3633: + x - the vector
3634: . m - first dimension of four dimensional array
3635: . n - second dimension of four dimensional array
3636: . p - third dimension of four dimensional array
3637: . q - fourth dimension of four dimensional array
3638: . mstart - first index you will use in first coordinate direction (often 0)
3639: . nstart - first index in the second coordinate direction (often 0)
3640: . pstart - first index in the third coordinate direction (often 0)
3641: - qstart - first index in the fourth coordinate direction (often 0)
3643: Output Parameter:
3644: . a - location to put pointer to the array
3646: Level: beginner
3648: Notes:
3649: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3650: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3651: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3652: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3654: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3656: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecPlaceArray()`,
3657: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3658: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3659: @*/
3660: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3661: {
3662: PetscInt i, N, j, k;
3663: const PetscScalar *aa;
3664: PetscScalar ***b, **c;
3666: PetscFunctionBegin;
3668: PetscAssertPointer(a, 10);
3670: PetscCall(VecGetLocalSize(x, &N));
3671: 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);
3672: PetscCall(VecGetArrayRead(x, &aa));
3674: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3675: b = (PetscScalar ***)((*a) + m);
3676: c = (PetscScalar **)(b + m * n);
3677: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3678: for (i = 0; i < m; i++)
3679: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3680: for (i = 0; i < m; i++)
3681: for (j = 0; j < n; j++)
3682: 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;
3683: *a -= mstart;
3684: PetscFunctionReturn(PETSC_SUCCESS);
3685: }
3687: /*@C
3688: VecRestoreArray4dRead - Restores a vector after `VecGetArray4d()` has been called.
3690: Logically Collective
3692: Input Parameters:
3693: + x - the vector
3694: . m - first dimension of four dimensional array
3695: . n - second dimension of the four dimensional array
3696: . p - third dimension of the four dimensional array
3697: . q - fourth dimension of the four dimensional array
3698: . mstart - first index you will use in first coordinate direction (often 0)
3699: . nstart - first index in the second coordinate direction (often 0)
3700: . pstart - first index in the third coordinate direction (often 0)
3701: . qstart - first index in the fourth coordinate direction (often 0)
3702: - a - location of pointer to array obtained from `VecGetArray4dRead()`
3704: Level: beginner
3706: Notes:
3707: For regular PETSc vectors this routine does not involve any copies. For
3708: any special vectors that do not store local vector data in a contiguous
3709: array, this routine will copy the data back into the underlying
3710: vector data structure from the array obtained with `VecGetArray()`.
3712: This routine actually zeros out the `a` pointer.
3714: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecPlaceArray()`,
3715: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`,
3716: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3717: @*/
3718: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3719: {
3720: void *dummy;
3722: PetscFunctionBegin;
3724: PetscAssertPointer(a, 10);
3726: dummy = (void *)(*a + mstart);
3727: PetscCall(PetscFree(dummy));
3728: PetscCall(VecRestoreArrayRead(x, NULL));
3729: *a = NULL;
3730: PetscFunctionReturn(PETSC_SUCCESS);
3731: }
3733: /*@
3734: VecLockGet - Get the current lock status of a vector
3736: Logically Collective
3738: Input Parameter:
3739: . x - the vector
3741: Output Parameter:
3742: . state - greater than zero indicates the vector is locked for read; less than zero indicates the vector is
3743: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3745: Level: advanced
3747: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3748: @*/
3749: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3750: {
3751: PetscFunctionBegin;
3753: PetscAssertPointer(state, 2);
3754: *state = x->lock;
3755: PetscFunctionReturn(PETSC_SUCCESS);
3756: }
3758: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3759: {
3760: PetscFunctionBegin;
3762: PetscAssertPointer(file, 2);
3763: PetscAssertPointer(func, 3);
3764: PetscAssertPointer(line, 4);
3765: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3766: {
3767: const int index = x->lockstack.currentsize - 1;
3769: *file = index < 0 ? NULL : x->lockstack.file[index];
3770: *func = index < 0 ? NULL : x->lockstack.function[index];
3771: *line = index < 0 ? 0 : x->lockstack.line[index];
3772: }
3773: #else
3774: *file = NULL;
3775: *func = NULL;
3776: *line = 0;
3777: #endif
3778: PetscFunctionReturn(PETSC_SUCCESS);
3779: }
3781: /*@
3782: VecLockReadPush - Push a read-only lock on a vector to prevent it from being written to
3784: Logically Collective
3786: Input Parameter:
3787: . x - the vector
3789: Level: intermediate
3791: Notes:
3792: If this is set then calls to `VecGetArray()` or `VecSetValues()` or any other routines that change the vectors values will generate an error.
3794: The call can be nested, i.e., called multiple times on the same vector, but each `VecLockReadPush()` has to have one matching
3795: `VecLockReadPop()`, which removes the latest read-only lock.
3797: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3798: @*/
3799: PetscErrorCode VecLockReadPush(Vec x)
3800: {
3801: PetscFunctionBegin;
3803: 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");
3804: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3805: {
3806: const char *file, *func;
3807: int index, line;
3809: if ((index = petscstack.currentsize - 2) < 0) {
3810: // vec was locked "outside" of petsc, either in user-land or main. the error message will
3811: // now show this function as the culprit, but it will include the stacktrace
3812: file = "unknown user-file";
3813: func = "unknown_user_function";
3814: line = 0;
3815: } else {
3816: file = petscstack.file[index];
3817: func = petscstack.function[index];
3818: line = petscstack.line[index];
3819: }
3820: PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
3821: }
3822: #endif
3823: PetscFunctionReturn(PETSC_SUCCESS);
3824: }
3826: /*@
3827: VecLockReadPop - Pop a read-only lock from a vector
3829: Logically Collective
3831: Input Parameter:
3832: . x - the vector
3834: Level: intermediate
3836: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
3837: @*/
3838: PetscErrorCode VecLockReadPop(Vec x)
3839: {
3840: PetscFunctionBegin;
3842: PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
3843: #if PetscDefined(USE_DEBUG) && !PetscDefined(HAVE_THREADSAFETY)
3844: {
3845: const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];
3847: PetscStackPop_Private(x->lockstack, previous);
3848: }
3849: #endif
3850: PetscFunctionReturn(PETSC_SUCCESS);
3851: }
3853: /*@
3854: VecLockWriteSet - Lock or unlock a vector for exclusive read/write access
3856: Logically Collective
3858: Input Parameters:
3859: + x - the vector
3860: - flg - `PETSC_TRUE` to lock the vector for exclusive read/write access; `PETSC_FALSE` to unlock it.
3862: Level: intermediate
3864: Notes:
3865: The function is useful in split-phase computations, which usually have a begin phase and an end phase.
3866: One can call `VecLockWriteSet`(x,`PETSC_TRUE`) in the begin phase to lock a vector for exclusive
3867: access, and call `VecLockWriteSet`(x,`PETSC_FALSE`) in the end phase to unlock the vector from exclusive
3868: access. In this way, one is ensured no other operations can access the vector in between. The code may like
3870: .vb
3871: VecGetArray(x,&xdata); // begin phase
3872: VecLockWriteSet(v,PETSC_TRUE);
3874: Other operations, which can not access x anymore (they can access xdata, of course)
3876: VecRestoreArray(x,&vdata); // end phase
3877: VecLockWriteSet(v,PETSC_FALSE);
3878: .ve
3880: The call can not be nested on the same vector, in other words, one can not call `VecLockWriteSet`(x,`PETSC_TRUE`)
3881: again before calling `VecLockWriteSet`(v,`PETSC_FALSE`).
3883: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
3884: @*/
3885: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
3886: {
3887: PetscFunctionBegin;
3889: if (flg) {
3890: 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");
3891: 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");
3892: x->lock = -1;
3893: } else {
3894: 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");
3895: x->lock = 0;
3896: }
3897: PetscFunctionReturn(PETSC_SUCCESS);
3898: }