Actual source code: rvector.c
1: /*
2: Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
3: These are the vector functions the user calls.
4: */
5: #include <petsc/private/vecimpl.h>
7: PetscInt VecGetSubVectorSavedStateId = -1;
9: #if PetscDefined(USE_DEBUG)
10: // this is a no-op '0' macro in optimized builds
11: PetscErrorCode VecValidValues_Internal(Vec vec, PetscInt argnum, PetscBool begin)
12: {
13: PetscFunctionBegin;
14: if (vec->petscnative || vec->ops->getarray) {
15: PetscInt n;
16: const PetscScalar *x;
17: PetscOffloadMask mask;
19: PetscCall(VecGetOffloadMask(vec, &mask));
20: if (!PetscOffloadHost(mask)) PetscFunctionReturn(PETSC_SUCCESS);
21: PetscCall(VecGetLocalSize(vec, &n));
22: PetscCall(VecGetArrayRead(vec, &x));
23: for (PetscInt i = 0; i < n; i++) {
24: if (begin) {
25: PetscCheck(!PetscIsInfOrNanScalar(x[i]), PETSC_COMM_SELF, PETSC_ERR_FP, "Vec entry at local location %" PetscInt_FMT " is not-a-number or infinite at beginning of function: Parameter number %" PetscInt_FMT, i, argnum);
26: } else {
27: PetscCheck(!PetscIsInfOrNanScalar(x[i]), PETSC_COMM_SELF, PETSC_ERR_FP, "Vec entry at local location %" PetscInt_FMT " is not-a-number or infinite at end of function: Parameter number %" PetscInt_FMT, i, argnum);
28: }
29: }
30: PetscCall(VecRestoreArrayRead(vec, &x));
31: }
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
34: #endif
36: /*@
37: VecMaxPointwiseDivide - Computes the maximum of the componentwise division `max = max_i abs(x[i]/y[i])`.
39: Logically Collective
41: Input Parameters:
42: + x - the numerators
43: - y - the denominators
45: Output Parameter:
46: . max - the result
48: Level: advanced
50: Notes:
51: `x` and `y` may be the same vector
53: if a particular `y[i]` is zero, it is treated as 1 in the above formula
55: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`
56: @*/
57: PetscErrorCode VecMaxPointwiseDivide(Vec x, Vec y, PetscReal *max)
58: {
59: PetscFunctionBegin;
62: PetscAssertPointer(max, 3);
65: PetscCheckSameTypeAndComm(x, 1, y, 2);
66: VecCheckSameSize(x, 1, y, 2);
67: VecCheckAssembled(x);
68: VecCheckAssembled(y);
69: PetscCall(VecLockReadPush(x));
70: PetscCall(VecLockReadPush(y));
71: PetscUseTypeMethod(x, maxpointwisedivide, y, max);
72: PetscCall(VecLockReadPop(x));
73: PetscCall(VecLockReadPop(y));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: /*@
78: VecDot - Computes the vector dot product.
80: Collective
82: Input Parameters:
83: + x - first vector
84: - y - second vector
86: Output Parameter:
87: . val - the dot product
89: Level: intermediate
91: Notes for Users of Complex Numbers:
92: For complex vectors, `VecDot()` computes
93: $ val = (x,y) = y^H x,
94: where y^H denotes the conjugate transpose of y. Note that this corresponds to the usual "mathematicians" complex
95: inner product where the SECOND argument gets the complex conjugate. Since the `BLASdot()` complex conjugates the first
96: first argument we call the `BLASdot()` with the arguments reversed.
98: Use `VecTDot()` for the indefinite form
99: $ val = (x,y) = y^T x,
100: where y^T denotes the transpose of y.
102: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDotRealPart()`
103: @*/
104: PetscErrorCode VecDot(Vec x, Vec y, PetscScalar *val)
105: {
106: PetscFunctionBegin;
109: PetscAssertPointer(val, 3);
112: PetscCheckSameTypeAndComm(x, 1, y, 2);
113: VecCheckSameSize(x, 1, y, 2);
114: VecCheckAssembled(x);
115: VecCheckAssembled(y);
117: PetscCall(VecLockReadPush(x));
118: PetscCall(VecLockReadPush(y));
119: PetscCall(PetscLogEventBegin(VEC_Dot, x, y, 0, 0));
120: PetscUseTypeMethod(x, dot, y, val);
121: PetscCall(PetscLogEventEnd(VEC_Dot, x, y, 0, 0));
122: PetscCall(VecLockReadPop(x));
123: PetscCall(VecLockReadPop(y));
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: /*@
128: VecDotRealPart - Computes the real part of the vector dot product.
130: Collective
132: Input Parameters:
133: + x - first vector
134: - y - second vector
136: Output Parameter:
137: . val - the real part of the dot product;
139: Level: intermediate
141: Notes for Users of Complex Numbers:
142: See `VecDot()` for more details on the definition of the dot product for complex numbers
144: For real numbers this returns the same value as `VecDot()`
146: For complex numbers in C^n (that is a vector of n components with a complex number for each component) this is equal to the usual real dot product on the
147: the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)
149: Developer Notes:
150: This is not currently optimized to compute only the real part of the dot product.
152: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDot()`, `VecDotNorm2()`
153: @*/
154: PetscErrorCode VecDotRealPart(Vec x, Vec y, PetscReal *val)
155: {
156: PetscScalar fdot;
158: PetscFunctionBegin;
159: PetscCall(VecDot(x, y, &fdot));
160: *val = PetscRealPart(fdot);
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: /*@
165: VecNorm - Computes the vector norm.
167: Collective
169: Input Parameters:
170: + x - the vector
171: - type - the type of the norm requested
173: Output Parameter:
174: . val - the norm
176: Level: intermediate
178: Notes:
179: See `NormType` for descriptions of each norm.
181: For complex numbers `NORM_1` will return the traditional 1 norm of the 2 norm of the complex
182: numbers; that is the 1 norm of the absolute values of the complex entries. In PETSc 3.6 and
183: earlier releases it returned the 1 norm of the 1 norm of the complex entries (what is
184: returned by the BLAS routine `asum()`). Both are valid norms but most people expect the former.
186: This routine stashes the computed norm value, repeated calls before the vector entries are
187: changed are then rapid since the precomputed value is immediately available. Certain vector
188: operations such as `VecSet()` store the norms so the value is immediately available and does
189: not need to be explicitly computed. `VecScale()` updates any stashed norm values, thus calls
190: after `VecScale()` do not need to explicitly recompute the norm.
192: .seealso: [](ch_vectors), `Vec`, `NormType`, `VecDot()`, `VecTDot()`, `VecDotBegin()`, `VecDotEnd()`, `VecNormAvailable()`,
193: `VecNormBegin()`, `VecNormEnd()`, `NormType()`
194: @*/
195: PetscErrorCode VecNorm(Vec x, NormType type, PetscReal *val)
196: {
197: PetscBool flg = PETSC_TRUE;
199: PetscFunctionBegin;
202: VecCheckAssembled(x);
204: PetscAssertPointer(val, 3);
206: PetscCall(VecNormAvailable(x, type, &flg, val));
207: // check that all MPI processes call this routine together and have same availability
208: if (PetscDefined(USE_DEBUG)) {
209: PetscMPIInt b0 = (PetscMPIInt)flg, b1[2], b2[2];
210: b1[0] = -b0;
211: b1[1] = b0;
212: PetscCall(MPIU_Allreduce(b1, b2, 2, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
213: PetscCheck(-b2[0] == b2[1], PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_WRONGSTATE, "Some MPI processes have cached %s norm, others do not. This may happen when some MPI processes call VecGetArray() and some others do not.", NormTypes[type]);
214: if (flg) {
215: PetscReal b1[2], b2[2];
216: b1[0] = -(*val);
217: b1[1] = *val;
218: PetscCall(MPIU_Allreduce(b1, b2, 2, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)x)));
219: PetscCheck(-b2[0] == b2[1], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Difference in cached %s norms: local %g", NormTypes[type], (double)*val);
220: }
221: }
222: if (flg) PetscFunctionReturn(PETSC_SUCCESS);
224: PetscCall(VecLockReadPush(x));
225: PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
226: PetscUseTypeMethod(x, norm, type, val);
227: PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
228: PetscCall(VecLockReadPop(x));
230: if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: /*@
235: VecNormAvailable - Returns the vector norm if it is already known. That is, it has been previously computed and cached in the vector
237: Not Collective
239: Input Parameters:
240: + x - the vector
241: - type - one of `NORM_1` (sum_i |x[i]|), `NORM_2` sqrt(sum_i (x[i])^2), `NORM_INFINITY` max_i |x[i]|. Also available
242: `NORM_1_AND_2`, which computes both norms and stores them
243: in a two element array.
245: Output Parameters:
246: + available - `PETSC_TRUE` if the val returned is valid
247: - val - the norm
249: Level: intermediate
251: Developer Notes:
252: `PETSC_HAVE_SLOW_BLAS_NORM2` will cause a C (loop unrolled) version of the norm to be used, rather
253: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
254: (as opposed to vendor provided) because the FORTRAN BLAS `NRM2()` routine is very slow.
256: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
257: `VecNormBegin()`, `VecNormEnd()`
258: @*/
259: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
260: {
261: PetscFunctionBegin;
264: PetscAssertPointer(available, 3);
265: PetscAssertPointer(val, 4);
267: if (type == NORM_1_AND_2) {
268: *available = PETSC_FALSE;
269: } else {
270: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
271: }
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: /*@
276: VecNormalize - Normalizes a vector by its 2-norm.
278: Collective
280: Input Parameter:
281: . x - the vector
283: Output Parameter:
284: . val - the vector norm before normalization. May be `NULL` if the value is not needed.
286: Level: intermediate
288: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
289: @*/
290: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
291: {
292: PetscReal norm;
294: PetscFunctionBegin;
297: PetscCall(VecSetErrorIfLocked(x, 1));
298: if (val) PetscAssertPointer(val, 2);
299: PetscCall(PetscLogEventBegin(VEC_Normalize, x, 0, 0, 0));
300: PetscCall(VecNorm(x, NORM_2, &norm));
301: if (norm == 0.0) PetscCall(PetscInfo(x, "Vector of zero norm can not be normalized; Returning only the zero norm\n"));
302: else if (PetscIsInfOrNanReal(norm)) PetscCall(PetscInfo(x, "Vector with Inf or Nan norm can not be normalized; Returning only the norm\n"));
303: else {
304: PetscScalar s = 1.0 / norm;
305: PetscCall(VecScale(x, s));
306: }
307: PetscCall(PetscLogEventEnd(VEC_Normalize, x, 0, 0, 0));
308: if (val) *val = norm;
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: /*@C
313: VecMax - Determines the vector component with maximum real part and its location.
315: Collective
317: Input Parameter:
318: . x - the vector
320: Output Parameters:
321: + p - the index of `val` (pass `NULL` if you don't want this) in the vector
322: - val - the maximum component
324: Level: intermediate
326: Notes:
327: Returns the value `PETSC_MIN_REAL` and negative `p` if the vector is of length 0.
329: Returns the smallest index with the maximum value
331: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
332: @*/
333: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
334: {
335: PetscFunctionBegin;
338: VecCheckAssembled(x);
339: if (p) PetscAssertPointer(p, 2);
340: PetscAssertPointer(val, 3);
341: PetscCall(VecLockReadPush(x));
342: PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
343: PetscUseTypeMethod(x, max, p, val);
344: PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
345: PetscCall(VecLockReadPop(x));
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: /*@C
350: VecMin - Determines the vector component with minimum real part and its location.
352: Collective
354: Input Parameter:
355: . x - the vector
357: Output Parameters:
358: + p - the index of `val` (pass `NULL` if you don't want this location) in the vector
359: - val - the minimum component
361: Level: intermediate
363: Notes:
364: Returns the value `PETSC_MAX_REAL` and negative `p` if the vector is of length 0.
366: This returns the smallest index with the minimum value
368: .seealso: [](ch_vectors), `Vec`, `VecMax()`
369: @*/
370: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
371: {
372: PetscFunctionBegin;
375: VecCheckAssembled(x);
376: if (p) PetscAssertPointer(p, 2);
377: PetscAssertPointer(val, 3);
378: PetscCall(VecLockReadPush(x));
379: PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
380: PetscUseTypeMethod(x, min, p, val);
381: PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
382: PetscCall(VecLockReadPop(x));
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: /*@
387: VecTDot - Computes an indefinite vector dot product. That is, this
388: routine does NOT use the complex conjugate.
390: Collective
392: Input Parameters:
393: + x - first vector
394: - y - second vector
396: Output Parameter:
397: . val - the dot product
399: Level: intermediate
401: Notes for Users of Complex Numbers:
402: For complex vectors, `VecTDot()` computes the indefinite form
403: $ val = (x,y) = y^T x,
404: where y^T denotes the transpose of y.
406: Use `VecDot()` for the inner product
407: $ val = (x,y) = y^H x,
408: where y^H denotes the conjugate transpose of y.
410: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
411: @*/
412: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
413: {
414: PetscFunctionBegin;
417: PetscAssertPointer(val, 3);
420: PetscCheckSameTypeAndComm(x, 1, y, 2);
421: VecCheckSameSize(x, 1, y, 2);
422: VecCheckAssembled(x);
423: VecCheckAssembled(y);
425: PetscCall(VecLockReadPush(x));
426: PetscCall(VecLockReadPush(y));
427: PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
428: PetscUseTypeMethod(x, tdot, y, val);
429: PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
430: PetscCall(VecLockReadPop(x));
431: PetscCall(VecLockReadPop(y));
432: PetscFunctionReturn(PETSC_SUCCESS);
433: }
435: PetscErrorCode VecScaleAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
436: {
437: PetscReal norms[4];
438: PetscBool flgs[4];
439: PetscScalar one = 1.0;
441: PetscFunctionBegin;
444: VecCheckAssembled(x);
445: PetscCall(VecSetErrorIfLocked(x, 1));
447: if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);
449: /* get current stashed norms */
450: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
452: PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
453: VecMethodDispatch(x, dctx, VecAsyncFnName(Scale), scale, (Vec, PetscScalar, PetscDeviceContext), alpha);
454: PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));
456: PetscCall(PetscObjectStateIncrease((PetscObject)x));
457: /* put the scaled stashed norms back into the Vec */
458: for (PetscInt i = 0; i < 4; i++) {
459: PetscReal ar = PetscAbsScalar(alpha);
460: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
461: }
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }
465: /*@
466: VecScale - Scales a vector.
468: Logically Collective
470: Input Parameters:
471: + x - the vector
472: - alpha - the scalar
474: Level: intermediate
476: Note:
477: For a vector with n components, `VecScale()` computes x[i] = alpha * x[i], for i=1,...,n.
479: .seealso: [](ch_vectors), `Vec`, `VecSet()`
480: @*/
481: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
482: {
483: PetscFunctionBegin;
484: PetscCall(VecScaleAsync_Private(x, alpha, NULL));
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: PetscErrorCode VecSetAsync_Private(Vec x, PetscScalar alpha, PetscDeviceContext dctx)
489: {
490: PetscFunctionBegin;
493: VecCheckAssembled(x);
495: PetscCall(VecSetErrorIfLocked(x, 1));
497: if (alpha == 0) {
498: PetscReal norm;
499: PetscBool set;
501: PetscCall(VecNormAvailable(x, NORM_2, &set, &norm));
502: if (set == PETSC_TRUE && norm == 0) PetscFunctionReturn(PETSC_SUCCESS);
503: }
504: PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
505: VecMethodDispatch(x, dctx, VecAsyncFnName(Set), set, (Vec, PetscScalar, PetscDeviceContext), alpha);
506: PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
507: PetscCall(PetscObjectStateIncrease((PetscObject)x));
509: /* norms can be simply set (if |alpha|*N not too large) */
510: {
511: PetscReal val = PetscAbsScalar(alpha);
512: const PetscInt N = x->map->N;
514: if (N == 0) {
515: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0));
516: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
517: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
518: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
519: } else if (val > PETSC_MAX_REAL / N) {
520: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
521: } else {
522: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
523: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
524: val *= PetscSqrtReal((PetscReal)N);
525: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
526: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
527: }
528: }
529: PetscFunctionReturn(PETSC_SUCCESS);
530: }
532: /*@
533: VecSet - Sets all components of a vector to a single scalar value.
535: Logically Collective
537: Input Parameters:
538: + x - the vector
539: - alpha - the scalar
541: Level: beginner
543: Notes:
544: For a vector of dimension n, `VecSet()` sets x[i] = alpha, for i=1,...,n,
545: so that all vector entries then equal the identical
546: scalar value, `alpha`. Use the more general routine
547: `VecSetValues()` to set different vector entries.
549: You CANNOT call this after you have called `VecSetValues()` but before you call
550: `VecAssemblyBegin()`
552: If `alpha` is zero and the norm of the vector is known to be zero then this skips the unneeded zeroing process
554: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
555: @*/
556: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
557: {
558: PetscFunctionBegin;
559: PetscCall(VecSetAsync_Private(x, alpha, NULL));
560: PetscFunctionReturn(PETSC_SUCCESS);
561: }
563: PetscErrorCode VecAXPYAsync_Private(Vec y, PetscScalar alpha, Vec x, PetscDeviceContext dctx)
564: {
565: PetscFunctionBegin;
570: PetscCheckSameTypeAndComm(x, 3, y, 1);
571: VecCheckSameSize(x, 3, y, 1);
572: VecCheckAssembled(x);
573: VecCheckAssembled(y);
575: if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
576: PetscCall(VecSetErrorIfLocked(y, 1));
577: if (x == y) {
578: PetscCall(VecScale(y, alpha + 1.0));
579: PetscFunctionReturn(PETSC_SUCCESS);
580: }
581: PetscCall(VecLockReadPush(x));
582: PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
583: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPY), axpy, (Vec, PetscScalar, Vec, PetscDeviceContext), alpha, x);
584: PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
585: PetscCall(VecLockReadPop(x));
586: PetscCall(PetscObjectStateIncrease((PetscObject)y));
587: PetscFunctionReturn(PETSC_SUCCESS);
588: }
589: /*@
590: VecAXPY - Computes `y = alpha x + y`.
592: Logically Collective
594: Input Parameters:
595: + alpha - the scalar
596: . x - vector scale by `alpha`
597: - y - vector accumulated into
599: Output Parameter:
600: . y - output vector
602: Level: intermediate
604: Notes:
605: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
606: .vb
607: VecAXPY(y,alpha,x) y = alpha x + y
608: VecAYPX(y,beta,x) y = x + beta y
609: VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
610: VecWAXPY(w,alpha,x,y) w = alpha x + y
611: VecAXPBYPCZ(z,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
612: VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
613: .ve
615: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
616: @*/
617: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
618: {
619: PetscFunctionBegin;
620: PetscCall(VecAXPYAsync_Private(y, alpha, x, NULL));
621: PetscFunctionReturn(PETSC_SUCCESS);
622: }
624: PetscErrorCode VecAYPXAsync_Private(Vec y, PetscScalar beta, Vec x, PetscDeviceContext dctx)
625: {
626: PetscFunctionBegin;
631: PetscCheckSameTypeAndComm(x, 3, y, 1);
632: VecCheckSameSize(x, 1, y, 3);
633: VecCheckAssembled(x);
634: VecCheckAssembled(y);
636: PetscCall(VecSetErrorIfLocked(y, 1));
637: if (x == y) {
638: PetscCall(VecScale(y, beta + 1.0));
639: PetscFunctionReturn(PETSC_SUCCESS);
640: }
641: PetscCall(VecLockReadPush(x));
642: if (beta == (PetscScalar)0.0) {
643: PetscCall(VecCopy(x, y));
644: } else {
645: PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
646: VecMethodDispatch(y, dctx, VecAsyncFnName(AYPX), aypx, (Vec, PetscScalar, Vec, PetscDeviceContext), beta, x);
647: PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
648: PetscCall(PetscObjectStateIncrease((PetscObject)y));
649: }
650: PetscCall(VecLockReadPop(x));
651: PetscFunctionReturn(PETSC_SUCCESS);
652: }
654: /*@
655: VecAYPX - Computes `y = x + beta y`.
657: Logically Collective
659: Input Parameters:
660: + beta - the scalar
661: . x - the unscaled vector
662: - y - the vector to be scaled
664: Output Parameter:
665: . y - output vector
667: Level: intermediate
669: Developer Notes:
670: The implementation is optimized for `beta` of -1.0, 0.0, and 1.0
672: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
673: @*/
674: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
675: {
676: PetscFunctionBegin;
677: PetscCall(VecAYPXAsync_Private(y, beta, x, NULL));
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: PetscErrorCode VecAXPBYAsync_Private(Vec y, PetscScalar alpha, PetscScalar beta, Vec x, PetscDeviceContext dctx)
682: {
683: PetscFunctionBegin;
688: PetscCheckSameTypeAndComm(x, 4, y, 1);
689: VecCheckSameSize(y, 1, x, 4);
690: VecCheckAssembled(x);
691: VecCheckAssembled(y);
694: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
695: if (x == y) {
696: PetscCall(VecScale(y, alpha + beta));
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: PetscCall(VecSetErrorIfLocked(y, 1));
701: PetscCall(VecLockReadPush(x));
702: PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
703: VecMethodDispatch(y, dctx, VecAsyncFnName(AXPBY), axpby, (Vec, PetscScalar, PetscScalar, Vec, PetscDeviceContext), alpha, beta, x);
704: PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
705: PetscCall(PetscObjectStateIncrease((PetscObject)y));
706: PetscCall(VecLockReadPop(x));
707: PetscFunctionReturn(PETSC_SUCCESS);
708: }
710: /*@
711: VecAXPBY - Computes `y = alpha x + beta y`.
713: Logically Collective
715: Input Parameters:
716: + alpha - first scalar
717: . beta - second scalar
718: . x - the first scaled vector
719: - y - the second scaled vector
721: Output Parameter:
722: . y - output vector
724: Level: intermediate
726: Developer Notes:
727: The implementation is optimized for `alpha` and/or `beta` values of 0.0 and 1.0
729: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
730: @*/
731: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
732: {
733: PetscFunctionBegin;
734: PetscCall(VecAXPBYAsync_Private(y, alpha, beta, x, NULL));
735: PetscFunctionReturn(PETSC_SUCCESS);
736: }
738: PetscErrorCode VecAXPBYPCZAsync_Private(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y, PetscDeviceContext dctx)
739: {
740: PetscFunctionBegin;
747: PetscCheckSameTypeAndComm(x, 5, y, 6);
748: PetscCheckSameTypeAndComm(x, 5, z, 1);
749: VecCheckSameSize(x, 5, y, 6);
750: VecCheckSameSize(x, 5, z, 1);
751: PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
752: PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
753: VecCheckAssembled(x);
754: VecCheckAssembled(y);
755: VecCheckAssembled(z);
759: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
761: PetscCall(VecSetErrorIfLocked(z, 1));
762: PetscCall(VecLockReadPush(x));
763: PetscCall(VecLockReadPush(y));
764: PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
765: VecMethodDispatch(z, dctx, VecAsyncFnName(AXPBYPCZ), axpbypcz, (Vec, PetscScalar, PetscScalar, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, beta, gamma, x, y);
766: PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
767: PetscCall(PetscObjectStateIncrease((PetscObject)z));
768: PetscCall(VecLockReadPop(x));
769: PetscCall(VecLockReadPop(y));
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
772: /*@
773: VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`
775: Logically Collective
777: Input Parameters:
778: + alpha - first scalar
779: . beta - second scalar
780: . gamma - third scalar
781: . x - first vector
782: . y - second vector
783: - z - third vector
785: Output Parameter:
786: . z - output vector
788: Level: intermediate
790: Note:
791: `x`, `y` and `z` must be different vectors
793: Developer Notes:
794: The implementation is optimized for `alpha` of 1.0 and `gamma` of 1.0 or 0.0
796: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
797: @*/
798: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
799: {
800: PetscFunctionBegin;
801: PetscCall(VecAXPBYPCZAsync_Private(z, alpha, beta, gamma, x, y, NULL));
802: PetscFunctionReturn(PETSC_SUCCESS);
803: }
805: PetscErrorCode VecWAXPYAsync_Private(Vec w, PetscScalar alpha, Vec x, Vec y, PetscDeviceContext dctx)
806: {
807: PetscFunctionBegin;
814: PetscCheckSameTypeAndComm(x, 3, y, 4);
815: PetscCheckSameTypeAndComm(y, 4, w, 1);
816: VecCheckSameSize(x, 3, y, 4);
817: VecCheckSameSize(x, 3, w, 1);
818: PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
819: PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
820: VecCheckAssembled(x);
821: VecCheckAssembled(y);
823: PetscCall(VecSetErrorIfLocked(w, 1));
825: PetscCall(VecLockReadPush(x));
826: PetscCall(VecLockReadPush(y));
827: if (alpha == (PetscScalar)0.0) {
828: PetscCall(VecCopyAsync_Private(y, w, dctx));
829: } else {
830: PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
831: VecMethodDispatch(w, dctx, VecAsyncFnName(WAXPY), waxpy, (Vec, PetscScalar, Vec, Vec, PetscDeviceContext), alpha, x, y);
832: PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
833: PetscCall(PetscObjectStateIncrease((PetscObject)w));
834: }
835: PetscCall(VecLockReadPop(x));
836: PetscCall(VecLockReadPop(y));
837: PetscFunctionReturn(PETSC_SUCCESS);
838: }
840: /*@
841: VecWAXPY - Computes `w = alpha x + y`.
843: Logically Collective
845: Input Parameters:
846: + alpha - the scalar
847: . x - first vector, multiplied by `alpha`
848: - y - second vector
850: Output Parameter:
851: . w - the result
853: Level: intermediate
855: Note:
856: `w` cannot be either `x` or `y`, but `x` and `y` can be the same
858: Developer Notes:
859: The implementation is optimized for alpha of -1.0, 0.0, and 1.0
861: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
862: @*/
863: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
864: {
865: PetscFunctionBegin;
866: PetscCall(VecWAXPYAsync_Private(w, alpha, x, y, NULL));
867: PetscFunctionReturn(PETSC_SUCCESS);
868: }
870: /*@C
871: VecSetValues - Inserts or adds values into certain locations of a vector.
873: Not Collective
875: Input Parameters:
876: + x - vector to insert in
877: . ni - number of elements to add
878: . ix - indices where to add
879: . y - array of values
880: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries
882: Level: beginner
884: Notes:
885: .vb
886: `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
887: .ve
889: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
890: options cannot be mixed without intervening calls to the assembly
891: routines.
893: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
894: MUST be called after all calls to `VecSetValues()` have been completed.
896: VecSetValues() uses 0-based indices in Fortran as well as in C.
898: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
899: negative indices may be passed in ix. These rows are
900: simply ignored. This allows easily inserting element load matrices
901: with homogeneous Dirichlet boundary conditions that you don't want represented
902: in the vector.
904: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
905: `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
906: @*/
907: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
908: {
909: PetscFunctionBeginHot;
911: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
912: PetscAssertPointer(ix, 3);
913: PetscAssertPointer(y, 4);
916: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
917: PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
918: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
919: PetscCall(PetscObjectStateIncrease((PetscObject)x));
920: PetscFunctionReturn(PETSC_SUCCESS);
921: }
923: /*@C
924: VecGetValues - Gets values from certain locations of a vector. Currently
925: can only get values on the same processor on which they are owned
927: Not Collective
929: Input Parameters:
930: + x - vector to get values from
931: . ni - number of elements to get
932: - ix - indices where to get them from (in global 1d numbering)
934: Output Parameter:
935: . y - array of values
937: Level: beginner
939: Notes:
940: The user provides the allocated array y; it is NOT allocated in this routine
942: `VecGetValues()` gets y[i] = x[ix[i]], for i=0,...,ni-1.
944: `VecAssemblyBegin()` and `VecAssemblyEnd()` MUST be called before calling this if `VecSetValues()` or related routine has been called
946: VecGetValues() uses 0-based indices in Fortran as well as in C.
948: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
949: negative indices may be passed in ix. These rows are
950: simply ignored.
952: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
953: @*/
954: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
955: {
956: PetscFunctionBegin;
958: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
959: PetscAssertPointer(ix, 3);
960: PetscAssertPointer(y, 4);
962: VecCheckAssembled(x);
963: PetscUseTypeMethod(x, getvalues, ni, ix, y);
964: PetscFunctionReturn(PETSC_SUCCESS);
965: }
967: /*@C
968: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
970: Not Collective
972: Input Parameters:
973: + x - vector to insert in
974: . ni - number of blocks to add
975: . ix - indices where to add in block count, rather than element count
976: . y - array of values
977: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries
979: Level: intermediate
981: Notes:
982: `VecSetValuesBlocked()` sets x[bs*ix[i]+j] = y[bs*i+j],
983: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
985: Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
986: options cannot be mixed without intervening calls to the assembly
987: routines.
989: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
990: MUST be called after all calls to `VecSetValuesBlocked()` have been completed.
992: `VecSetValuesBlocked()` uses 0-based indices in Fortran as well as in C.
994: Negative indices may be passed in ix, these rows are
995: simply ignored. This allows easily inserting element load matrices
996: with homogeneous Dirichlet boundary conditions that you don't want represented
997: in the vector.
999: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
1000: `VecSetValues()`
1001: @*/
1002: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1003: {
1004: PetscFunctionBeginHot;
1006: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1007: PetscAssertPointer(ix, 3);
1008: PetscAssertPointer(y, 4);
1011: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1012: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1013: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1014: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1015: PetscFunctionReturn(PETSC_SUCCESS);
1016: }
1018: /*@C
1019: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
1020: using a local ordering of the nodes.
1022: Not Collective
1024: Input Parameters:
1025: + x - vector to insert in
1026: . ni - number of elements to add
1027: . ix - indices where to add
1028: . y - array of values
1029: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1031: Level: intermediate
1033: Notes:
1034: `VecSetValuesLocal()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
1036: Calls to `VecSetValuesLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1037: options cannot be mixed without intervening calls to the assembly
1038: routines.
1040: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1041: MUST be called after all calls to `VecSetValuesLocal()` have been completed.
1043: `VecSetValuesLocal()` uses 0-based indices in Fortran as well as in C.
1045: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1046: `VecSetValuesBlockedLocal()`
1047: @*/
1048: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1049: {
1050: PetscInt lixp[128], *lix = lixp;
1052: PetscFunctionBeginHot;
1054: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1055: PetscAssertPointer(ix, 3);
1056: PetscAssertPointer(y, 4);
1059: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1060: if (!x->ops->setvalueslocal) {
1061: if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
1062: if (x->map->mapping) {
1063: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1064: PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1065: PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1066: if (ni > 128) PetscCall(PetscFree(lix));
1067: } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1068: } else PetscUseTypeMethod(x, setvalueslocal, ni, ix, y, iora);
1069: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1070: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1071: PetscFunctionReturn(PETSC_SUCCESS);
1072: }
1074: /*@
1075: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1076: using a local ordering of the nodes.
1078: Not Collective
1080: Input Parameters:
1081: + x - vector to insert in
1082: . ni - number of blocks to add
1083: . ix - indices where to add in block count, not element count
1084: . y - array of values
1085: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1087: Level: intermediate
1089: Notes:
1090: `VecSetValuesBlockedLocal()` sets x[bs*ix[i]+j] = y[bs*i+j],
1091: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with `VecSetBlockSize()`.
1093: Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1094: options cannot be mixed without intervening calls to the assembly
1095: routines.
1097: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1098: MUST be called after all calls to `VecSetValuesBlockedLocal()` have been completed.
1100: `VecSetValuesBlockedLocal()` uses 0-based indices in Fortran as well as in C.
1102: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1103: `VecSetLocalToGlobalMapping()`
1104: @*/
1105: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1106: {
1107: PetscInt lixp[128], *lix = lixp;
1109: PetscFunctionBeginHot;
1111: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1112: PetscAssertPointer(ix, 3);
1113: PetscAssertPointer(y, 4);
1115: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1116: if (PetscUnlikely(!x->map->mapping && x->ops->getlocaltoglobalmapping)) PetscUseTypeMethod(x, getlocaltoglobalmapping, &x->map->mapping);
1117: if (x->map->mapping) {
1118: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1119: PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1120: PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1121: if (ni > 128) PetscCall(PetscFree(lix));
1122: } else {
1123: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1124: }
1125: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1126: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1127: PetscFunctionReturn(PETSC_SUCCESS);
1128: }
1130: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1131: {
1132: PetscFunctionBegin;
1135: VecCheckAssembled(x);
1137: if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1138: PetscAssertPointer(y, 3);
1139: for (PetscInt i = 0; i < nv; ++i) {
1142: PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1143: VecCheckSameSize(x, 1, y[i], 3);
1144: VecCheckAssembled(y[i]);
1145: PetscCall(VecLockReadPush(y[i]));
1146: }
1147: PetscAssertPointer(result, 4);
1150: PetscCall(VecLockReadPush(x));
1151: PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1152: PetscCall((*mxdot)(x, nv, y, result));
1153: PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1154: PetscCall(VecLockReadPop(x));
1155: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1156: PetscFunctionReturn(PETSC_SUCCESS);
1157: }
1159: /*@
1160: VecMTDot - Computes indefinite vector multiple dot products.
1161: That is, it does NOT use the complex conjugate.
1163: Collective
1165: Input Parameters:
1166: + x - one vector
1167: . nv - number of vectors
1168: - y - array of vectors. Note that vectors are pointers
1170: Output Parameter:
1171: . val - array of the dot products
1173: Level: intermediate
1175: Notes for Users of Complex Numbers:
1176: For complex vectors, `VecMTDot()` computes the indefinite form
1177: $ val = (x,y) = y^T x,
1178: where y^T denotes the transpose of y.
1180: Use `VecMDot()` for the inner product
1181: $ val = (x,y) = y^H x,
1182: where y^H denotes the conjugate transpose of y.
1184: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1185: @*/
1186: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1187: {
1188: PetscFunctionBegin;
1190: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1191: PetscFunctionReturn(PETSC_SUCCESS);
1192: }
1194: /*@
1195: VecMDot - Computes multiple vector dot products.
1197: Collective
1199: Input Parameters:
1200: + x - one vector
1201: . nv - number of vectors
1202: - y - array of vectors.
1204: Output Parameter:
1205: . val - array of the dot products (does not allocate the array)
1207: Level: intermediate
1209: Notes for Users of Complex Numbers:
1210: For complex vectors, `VecMDot()` computes
1211: $ val = (x,y) = y^H x,
1212: where y^H denotes the conjugate transpose of y.
1214: Use `VecMTDot()` for the indefinite form
1215: $ val = (x,y) = y^T x,
1216: where y^T denotes the transpose of y.
1218: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`
1219: @*/
1220: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1221: {
1222: PetscFunctionBegin;
1224: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1225: PetscFunctionReturn(PETSC_SUCCESS);
1226: }
1228: PetscErrorCode VecMAXPYAsync_Private(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[], PetscDeviceContext dctx)
1229: {
1230: PetscFunctionBegin;
1232: VecCheckAssembled(y);
1234: PetscCall(VecSetErrorIfLocked(y, 1));
1235: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1236: if (nv) {
1237: PetscInt zeros = 0;
1239: PetscAssertPointer(alpha, 3);
1240: PetscAssertPointer(x, 4);
1241: for (PetscInt i = 0; i < nv; ++i) {
1245: PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1246: VecCheckSameSize(y, 1, x[i], 4);
1247: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1248: VecCheckAssembled(x[i]);
1249: PetscCall(VecLockReadPush(x[i]));
1250: zeros += alpha[i] == (PetscScalar)0.0;
1251: }
1253: if (zeros < nv) {
1254: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1255: VecMethodDispatch(y, dctx, VecAsyncFnName(MAXPY), maxpy, (Vec, PetscInt, const PetscScalar[], Vec[], PetscDeviceContext), nv, alpha, x);
1256: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1257: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1258: }
1260: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1261: }
1262: PetscFunctionReturn(PETSC_SUCCESS);
1263: }
1265: /*@
1266: VecMAXPY - Computes `y = y + sum alpha[i] x[i]`
1268: Logically Collective
1270: Input Parameters:
1271: + nv - number of scalars and x-vectors
1272: . alpha - array of scalars
1273: . y - one vector
1274: - x - array of vectors
1276: Level: intermediate
1278: Note:
1279: `y` cannot be any of the `x` vectors
1281: .seealso: [](ch_vectors), `Vec`, `VecMAXPBY()`,`VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1282: @*/
1283: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1284: {
1285: PetscFunctionBegin;
1286: PetscCall(VecMAXPYAsync_Private(y, nv, alpha, x, NULL));
1287: PetscFunctionReturn(PETSC_SUCCESS);
1288: }
1290: /*@
1291: VecMAXPBY - Computes `y = beta y + sum alpha[i] x[i]`
1293: Logically Collective
1295: Input Parameters:
1296: + nv - number of scalars and x-vectors
1297: . alpha - array of scalars
1298: . beta - scalar
1299: . y - one vector
1300: - x - array of vectors
1302: Level: intermediate
1304: Note:
1305: `y` cannot be any of the `x` vectors.
1307: Developer Notes:
1308: This is a convenience routine, but implementations might be able to optimize it, for example, when `beta` is zero.
1310: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1311: @*/
1312: PetscErrorCode VecMAXPBY(Vec y, PetscInt nv, const PetscScalar alpha[], PetscScalar beta, Vec x[])
1313: {
1314: PetscFunctionBegin;
1316: VecCheckAssembled(y);
1318: PetscCall(VecSetErrorIfLocked(y, 1));
1319: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1322: if (y->ops->maxpby) {
1323: PetscInt zeros = 0;
1325: if (nv) {
1326: PetscAssertPointer(alpha, 3);
1327: PetscAssertPointer(x, 5);
1328: }
1330: for (PetscInt i = 0; i < nv; ++i) { // scan all alpha[]
1334: PetscCheckSameTypeAndComm(y, 1, x[i], 5);
1335: VecCheckSameSize(y, 1, x[i], 5);
1336: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1337: VecCheckAssembled(x[i]);
1338: PetscCall(VecLockReadPush(x[i]));
1339: zeros += alpha[i] == (PetscScalar)0.0;
1340: }
1342: if (zeros < nv) { // has nonzero alpha
1343: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1344: PetscUseTypeMethod(y, maxpby, nv, alpha, beta, x);
1345: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1346: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1347: } else {
1348: PetscCall(VecScale(y, beta));
1349: }
1351: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1352: } else { // no maxpby
1353: if (beta == 0.0) PetscCall(VecSet(y, 0.0));
1354: else PetscCall(VecScale(y, beta));
1355: PetscCall(VecMAXPY(y, nv, alpha, x));
1356: }
1357: PetscFunctionReturn(PETSC_SUCCESS);
1358: }
1360: /*@
1361: VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
1362: in the order they appear in the array. The concatenated vector resides on the same
1363: communicator and is the same type as the source vectors.
1365: Collective
1367: Input Parameters:
1368: + nx - number of vectors to be concatenated
1369: - X - array containing the vectors to be concatenated in the order of concatenation
1371: Output Parameters:
1372: + Y - concatenated vector
1373: - x_is - array of index sets corresponding to the concatenated components of `Y` (pass `NULL` if not needed)
1375: Level: advanced
1377: Notes:
1378: Concatenation is similar to the functionality of a `VECNEST` object; they both represent combination of
1379: different vector spaces. However, concatenated vectors do not store any information about their
1380: sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1381: manipulation of data in the concatenated vector that corresponds to the original components at creation.
1383: This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1384: has to operate on combined vector spaces and cannot utilize `VECNEST` objects due to incompatibility with
1385: bound projections.
1387: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1388: @*/
1389: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1390: {
1391: MPI_Comm comm;
1392: VecType vec_type;
1393: Vec Ytmp, Xtmp;
1394: IS *is_tmp;
1395: PetscInt i, shift = 0, Xnl, Xng, Xbegin;
1397: PetscFunctionBegin;
1401: PetscAssertPointer(Y, 3);
1403: if ((*X)->ops->concatenate) {
1404: /* use the dedicated concatenation function if available */
1405: PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1406: } else {
1407: /* loop over vectors and start creating IS */
1408: comm = PetscObjectComm((PetscObject)*X);
1409: PetscCall(VecGetType(*X, &vec_type));
1410: PetscCall(PetscMalloc1(nx, &is_tmp));
1411: for (i = 0; i < nx; i++) {
1412: PetscCall(VecGetSize(X[i], &Xng));
1413: PetscCall(VecGetLocalSize(X[i], &Xnl));
1414: PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1415: PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1416: shift += Xng;
1417: }
1418: /* create the concatenated vector */
1419: PetscCall(VecCreate(comm, &Ytmp));
1420: PetscCall(VecSetType(Ytmp, vec_type));
1421: PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1422: PetscCall(VecSetUp(Ytmp));
1423: /* copy data from X array to Y and return */
1424: for (i = 0; i < nx; i++) {
1425: PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1426: PetscCall(VecCopy(X[i], Xtmp));
1427: PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1428: }
1429: *Y = Ytmp;
1430: if (x_is) {
1431: *x_is = is_tmp;
1432: } else {
1433: for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1434: PetscCall(PetscFree(is_tmp));
1435: }
1436: }
1437: PetscFunctionReturn(PETSC_SUCCESS);
1438: }
1440: /* A helper function for VecGetSubVector to check if we can implement it with no-copy (i.e. the subvector shares
1441: memory with the original vector), and the block size of the subvector.
1443: Input Parameters:
1444: + X - the original vector
1445: - is - the index set of the subvector
1447: Output Parameters:
1448: + contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
1449: . start - start of contiguous block, as an offset from the start of the ownership range of the original vector
1450: - blocksize - the block size of the subvector
1452: */
1453: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1454: {
1455: PetscInt gstart, gend, lstart;
1456: PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1457: PetscInt n, N, ibs, vbs, bs = -1;
1459: PetscFunctionBegin;
1460: PetscCall(ISGetLocalSize(is, &n));
1461: PetscCall(ISGetSize(is, &N));
1462: PetscCall(ISGetBlockSize(is, &ibs));
1463: PetscCall(VecGetBlockSize(X, &vbs));
1464: PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1465: PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1466: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1467: if (ibs > 1) {
1468: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1469: bs = ibs;
1470: } else {
1471: if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1472: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1473: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1474: }
1476: *contig = red[0];
1477: *start = lstart;
1478: *blocksize = bs;
1479: PetscFunctionReturn(PETSC_SUCCESS);
1480: }
1482: /* A helper function for VecGetSubVector, to be used when we have to build a standalone subvector through VecScatter
1484: Input Parameters:
1485: + X - the original vector
1486: . is - the index set of the subvector
1487: - bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()
1489: Output Parameter:
1490: . Z - the subvector, which will compose the VecScatter context on output
1491: */
1492: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1493: {
1494: PetscInt n, N;
1495: VecScatter vscat;
1496: Vec Y;
1498: PetscFunctionBegin;
1499: PetscCall(ISGetLocalSize(is, &n));
1500: PetscCall(ISGetSize(is, &N));
1501: PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1502: PetscCall(VecSetSizes(Y, n, N));
1503: PetscCall(VecSetBlockSize(Y, bs));
1504: PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1505: PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1506: PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1507: PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1508: PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1509: PetscCall(VecScatterDestroy(&vscat));
1510: *Z = Y;
1511: PetscFunctionReturn(PETSC_SUCCESS);
1512: }
1514: /*@
1515: VecGetSubVector - Gets a vector representing part of another vector
1517: Collective
1519: Input Parameters:
1520: + X - vector from which to extract a subvector
1521: - is - index set representing portion of `X` to extract
1523: Output Parameter:
1524: . Y - subvector corresponding to `is`
1526: Level: advanced
1528: Notes:
1529: The subvector `Y` should be returned with `VecRestoreSubVector()`.
1530: `X` and `is` must be defined on the same communicator
1532: Changes to the subvector will be reflected in the `X` vector on the call to `VecRestoreSubVector()`.
1534: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1535: modifying the subvector. Other non-overlapping subvectors can still be obtained from `X` using this function.
1537: 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`.
1539: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1540: @*/
1541: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1542: {
1543: Vec Z;
1545: PetscFunctionBegin;
1548: PetscCheckSameComm(X, 1, is, 2);
1549: PetscAssertPointer(Y, 3);
1550: if (X->ops->getsubvector) {
1551: PetscUseTypeMethod(X, getsubvector, is, &Z);
1552: } else { /* Default implementation currently does no caching */
1553: PetscBool contig;
1554: PetscInt n, N, start, bs;
1556: PetscCall(ISGetLocalSize(is, &n));
1557: PetscCall(ISGetSize(is, &N));
1558: PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1559: if (contig) { /* We can do a no-copy implementation */
1560: const PetscScalar *x;
1561: PetscInt state = 0;
1562: PetscBool isstd, iscuda, iship;
1564: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1565: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1566: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1567: if (iscuda) {
1568: #if defined(PETSC_HAVE_CUDA)
1569: const PetscScalar *x_d;
1570: PetscMPIInt size;
1571: PetscOffloadMask flg;
1573: PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1574: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1575: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1576: if (x) x += start;
1577: if (x_d) x_d += start;
1578: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1579: if (size == 1) {
1580: PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1581: } else {
1582: PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1583: }
1584: Z->offloadmask = flg;
1585: #endif
1586: } else if (iship) {
1587: #if defined(PETSC_HAVE_HIP)
1588: const PetscScalar *x_d;
1589: PetscMPIInt size;
1590: PetscOffloadMask flg;
1592: PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1593: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1594: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1595: if (x) x += start;
1596: if (x_d) x_d += start;
1597: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1598: if (size == 1) {
1599: PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1600: } else {
1601: PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1602: }
1603: Z->offloadmask = flg;
1604: #endif
1605: } else if (isstd) {
1606: PetscMPIInt size;
1608: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1609: PetscCall(VecGetArrayRead(X, &x));
1610: if (x) x += start;
1611: if (size == 1) {
1612: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1613: } else {
1614: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1615: }
1616: PetscCall(VecRestoreArrayRead(X, &x));
1617: } else { /* default implementation: use place array */
1618: PetscCall(VecGetArrayRead(X, &x));
1619: PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1620: PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1621: PetscCall(VecSetSizes(Z, n, N));
1622: PetscCall(VecSetBlockSize(Z, bs));
1623: PetscCall(VecPlaceArray(Z, PetscSafePointerPlusOffset(x, start)));
1624: PetscCall(VecRestoreArrayRead(X, &x));
1625: }
1627: /* this is relevant only in debug mode */
1628: PetscCall(VecLockGet(X, &state));
1629: if (state) PetscCall(VecLockReadPush(Z));
1630: Z->ops->placearray = NULL;
1631: Z->ops->replacearray = NULL;
1632: } else { /* Have to create a scatter and do a copy */
1633: PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1634: }
1635: }
1636: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1637: if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1638: PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1639: *Y = Z;
1640: PetscFunctionReturn(PETSC_SUCCESS);
1641: }
1643: /*@
1644: VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`
1646: Collective
1648: Input Parameters:
1649: + X - vector from which subvector was obtained
1650: . is - index set representing the subset of `X`
1651: - Y - subvector being restored
1653: Level: advanced
1655: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1656: @*/
1657: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1658: {
1659: PETSC_UNUSED PetscObjectState dummystate = 0;
1660: PetscBool unchanged;
1662: PetscFunctionBegin;
1665: PetscCheckSameComm(X, 1, is, 2);
1666: PetscAssertPointer(Y, 3);
1669: if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1670: else {
1671: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1672: if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1673: VecScatter scatter;
1674: PetscInt state;
1676: PetscCall(VecLockGet(X, &state));
1677: PetscCheck(state == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec X is locked for read-only or read/write access");
1679: PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1680: if (scatter) {
1681: PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1682: PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1683: } else {
1684: PetscBool iscuda, iship;
1685: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1686: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1688: if (iscuda) {
1689: #if defined(PETSC_HAVE_CUDA)
1690: PetscOffloadMask ymask = (*Y)->offloadmask;
1692: /* The offloadmask of X dictates where to move memory
1693: If X GPU data is valid, then move Y data on GPU if needed
1694: Otherwise, move back to the CPU */
1695: switch (X->offloadmask) {
1696: case PETSC_OFFLOAD_BOTH:
1697: if (ymask == PETSC_OFFLOAD_CPU) {
1698: PetscCall(VecCUDAResetArray(*Y));
1699: } else if (ymask == PETSC_OFFLOAD_GPU) {
1700: X->offloadmask = PETSC_OFFLOAD_GPU;
1701: }
1702: break;
1703: case PETSC_OFFLOAD_GPU:
1704: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1705: break;
1706: case PETSC_OFFLOAD_CPU:
1707: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1708: break;
1709: case PETSC_OFFLOAD_UNALLOCATED:
1710: case PETSC_OFFLOAD_KOKKOS:
1711: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1712: }
1713: #endif
1714: } else if (iship) {
1715: #if defined(PETSC_HAVE_HIP)
1716: PetscOffloadMask ymask = (*Y)->offloadmask;
1718: /* The offloadmask of X dictates where to move memory
1719: If X GPU data is valid, then move Y data on GPU if needed
1720: Otherwise, move back to the CPU */
1721: switch (X->offloadmask) {
1722: case PETSC_OFFLOAD_BOTH:
1723: if (ymask == PETSC_OFFLOAD_CPU) {
1724: PetscCall(VecHIPResetArray(*Y));
1725: } else if (ymask == PETSC_OFFLOAD_GPU) {
1726: X->offloadmask = PETSC_OFFLOAD_GPU;
1727: }
1728: break;
1729: case PETSC_OFFLOAD_GPU:
1730: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1731: break;
1732: case PETSC_OFFLOAD_CPU:
1733: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1734: break;
1735: case PETSC_OFFLOAD_UNALLOCATED:
1736: case PETSC_OFFLOAD_KOKKOS:
1737: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1738: }
1739: #endif
1740: } else {
1741: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1742: PetscCall(VecResetArray(*Y));
1743: }
1744: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1745: }
1746: }
1747: }
1748: PetscCall(VecDestroy(Y));
1749: PetscFunctionReturn(PETSC_SUCCESS);
1750: }
1752: /*@
1753: VecCreateLocalVector - Creates a vector object suitable for use with `VecGetLocalVector()` and friends. You must call `VecDestroy()` when the
1754: vector is no longer needed.
1756: Not Collective.
1758: Input Parameter:
1759: . v - The vector for which the local vector is desired.
1761: Output Parameter:
1762: . w - Upon exit this contains the local vector.
1764: Level: beginner
1766: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1767: @*/
1768: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1769: {
1770: PetscMPIInt size;
1772: PetscFunctionBegin;
1774: PetscAssertPointer(w, 2);
1775: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1776: if (size == 1) PetscCall(VecDuplicate(v, w));
1777: else if (v->ops->createlocalvector) PetscUseTypeMethod(v, createlocalvector, w);
1778: else {
1779: VecType type;
1780: PetscInt n;
1782: PetscCall(VecCreate(PETSC_COMM_SELF, w));
1783: PetscCall(VecGetLocalSize(v, &n));
1784: PetscCall(VecSetSizes(*w, n, n));
1785: PetscCall(VecGetBlockSize(v, &n));
1786: PetscCall(VecSetBlockSize(*w, n));
1787: PetscCall(VecGetType(v, &type));
1788: PetscCall(VecSetType(*w, type));
1789: }
1790: PetscFunctionReturn(PETSC_SUCCESS);
1791: }
1793: /*@
1794: VecGetLocalVectorRead - Maps the local portion of a vector into a
1795: vector.
1797: Not Collective.
1799: Input Parameter:
1800: . v - The vector for which the local vector is desired.
1802: Output Parameter:
1803: . w - Upon exit this contains the local vector.
1805: Level: beginner
1807: Notes:
1808: You must call `VecRestoreLocalVectorRead()` when the local
1809: vector is no longer needed.
1811: This function is similar to `VecGetArrayRead()` which maps the local
1812: portion into a raw pointer. `VecGetLocalVectorRead()` is usually
1813: almost as efficient as `VecGetArrayRead()` but in certain circumstances
1814: `VecGetLocalVectorRead()` can be much more efficient than
1815: `VecGetArrayRead()`. This is because the construction of a contiguous
1816: array representing the vector data required by `VecGetArrayRead()` can
1817: be an expensive operation for certain vector types. For example, for
1818: GPU vectors `VecGetArrayRead()` requires that the data between device
1819: and host is synchronized.
1821: Unlike `VecGetLocalVector()`, this routine is not collective and
1822: preserves cached information.
1824: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1825: @*/
1826: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1827: {
1828: PetscFunctionBegin;
1831: VecCheckSameLocalSize(v, 1, w, 2);
1832: if (v->ops->getlocalvectorread) {
1833: PetscUseTypeMethod(v, getlocalvectorread, w);
1834: } else {
1835: PetscScalar *a;
1837: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1838: PetscCall(VecPlaceArray(w, a));
1839: }
1840: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1841: PetscCall(VecLockReadPush(v));
1842: PetscCall(VecLockReadPush(w));
1843: PetscFunctionReturn(PETSC_SUCCESS);
1844: }
1846: /*@
1847: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1848: previously mapped into a vector using `VecGetLocalVectorRead()`.
1850: Not Collective.
1852: Input Parameters:
1853: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
1854: - w - The vector into which the local portion of `v` was mapped.
1856: Level: beginner
1858: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1859: @*/
1860: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1861: {
1862: PetscFunctionBegin;
1865: if (v->ops->restorelocalvectorread) {
1866: PetscUseTypeMethod(v, restorelocalvectorread, w);
1867: } else {
1868: const PetscScalar *a;
1870: PetscCall(VecGetArrayRead(w, &a));
1871: PetscCall(VecRestoreArrayRead(v, &a));
1872: PetscCall(VecResetArray(w));
1873: }
1874: PetscCall(VecLockReadPop(v));
1875: PetscCall(VecLockReadPop(w));
1876: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1877: PetscFunctionReturn(PETSC_SUCCESS);
1878: }
1880: /*@
1881: VecGetLocalVector - Maps the local portion of a vector into a
1882: vector.
1884: Collective
1886: Input Parameter:
1887: . v - The vector for which the local vector is desired.
1889: Output Parameter:
1890: . w - Upon exit this contains the local vector.
1892: Level: beginner
1894: Notes:
1895: You must call `VecRestoreLocalVector()` when the local
1896: vector is no longer needed.
1898: This function is similar to `VecGetArray()` which maps the local
1899: portion into a raw pointer. `VecGetLocalVector()` is usually about as
1900: efficient as `VecGetArray()` but in certain circumstances
1901: `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1902: This is because the construction of a contiguous array representing
1903: the vector data required by `VecGetArray()` can be an expensive
1904: operation for certain vector types. For example, for GPU vectors
1905: `VecGetArray()` requires that the data between device and host is
1906: synchronized.
1908: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1909: @*/
1910: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1911: {
1912: PetscFunctionBegin;
1915: VecCheckSameLocalSize(v, 1, w, 2);
1916: if (v->ops->getlocalvector) {
1917: PetscUseTypeMethod(v, getlocalvector, w);
1918: } else {
1919: PetscScalar *a;
1921: PetscCall(VecGetArray(v, &a));
1922: PetscCall(VecPlaceArray(w, a));
1923: }
1924: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1925: PetscFunctionReturn(PETSC_SUCCESS);
1926: }
1928: /*@
1929: VecRestoreLocalVector - Unmaps the local portion of a vector
1930: previously mapped into a vector using `VecGetLocalVector()`.
1932: Logically Collective.
1934: Input Parameters:
1935: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
1936: - w - The vector into which the local portion of `v` was mapped.
1938: Level: beginner
1940: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1941: @*/
1942: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1943: {
1944: PetscFunctionBegin;
1947: if (v->ops->restorelocalvector) {
1948: PetscUseTypeMethod(v, restorelocalvector, w);
1949: } else {
1950: PetscScalar *a;
1951: PetscCall(VecGetArray(w, &a));
1952: PetscCall(VecRestoreArray(v, &a));
1953: PetscCall(VecResetArray(w));
1954: }
1955: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1956: PetscCall(PetscObjectStateIncrease((PetscObject)v));
1957: PetscFunctionReturn(PETSC_SUCCESS);
1958: }
1960: /*@C
1961: VecGetArray - Returns a pointer to a contiguous array that contains this
1962: MPI processes's portion of the vector data
1964: Logically Collective
1966: Input Parameter:
1967: . x - the vector
1969: Output Parameter:
1970: . a - location to put pointer to the array
1972: Level: beginner
1974: Notes:
1975: For the standard PETSc vectors, `VecGetArray()` returns a pointer to the local data array and
1976: does not use any copies. If the underlying vector data is not stored in a contiguous array
1977: this routine will copy the data to a contiguous array and return a pointer to that. You MUST
1978: call `VecRestoreArray()` when you no longer need access to the array.
1980: Fortran Notes:
1981: `VecGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayF90()`
1983: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
1984: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
1985: @*/
1986: PetscErrorCode VecGetArray(Vec x, PetscScalar **a)
1987: {
1988: PetscFunctionBegin;
1990: PetscCall(VecSetErrorIfLocked(x, 1));
1991: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
1992: PetscUseTypeMethod(x, getarray, a);
1993: } else if (x->petscnative) { /* VECSTANDARD */
1994: *a = *((PetscScalar **)x->data);
1995: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
1996: PetscFunctionReturn(PETSC_SUCCESS);
1997: }
1999: /*@C
2000: VecRestoreArray - Restores a vector after `VecGetArray()` has been called and the array is no longer needed
2002: Logically Collective
2004: Input Parameters:
2005: + x - the vector
2006: - a - location of pointer to array obtained from `VecGetArray()`
2008: Level: beginner
2010: Fortran Notes:
2011: `VecRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayF90()`
2013: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2014: `VecGetArrayPair()`, `VecRestoreArrayPair()`
2015: @*/
2016: PetscErrorCode VecRestoreArray(Vec x, PetscScalar **a)
2017: {
2018: PetscFunctionBegin;
2020: if (a) PetscAssertPointer(a, 2);
2021: if (x->ops->restorearray) {
2022: PetscUseTypeMethod(x, restorearray, a);
2023: } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
2024: if (a) *a = NULL;
2025: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2026: PetscFunctionReturn(PETSC_SUCCESS);
2027: }
2028: /*@C
2029: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
2031: Not Collective
2033: Input Parameter:
2034: . x - the vector
2036: Output Parameter:
2037: . a - the array
2039: Level: beginner
2041: Notes:
2042: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2044: Unlike `VecGetArray()`, preserves cached information like vector norms.
2046: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
2047: implementations may require a copy, but such implementations should cache the contiguous representation so that
2048: only one copy is performed when this routine is called multiple times in sequence.
2050: Fortran Notes:
2051: `VecGetArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayReadF90()`
2053: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2054: @*/
2055: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar **a)
2056: {
2057: PetscFunctionBegin;
2059: PetscAssertPointer(a, 2);
2060: if (x->ops->getarrayread) {
2061: PetscUseTypeMethod(x, getarrayread, a);
2062: } else if (x->ops->getarray) {
2063: PetscObjectState state;
2065: /* VECNEST, VECCUDA, VECKOKKOS etc */
2066: // x->ops->getarray may bump the object state, but since we know this is a read-only get
2067: // we can just undo that
2068: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2069: PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
2070: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2071: } else if (x->petscnative) {
2072: /* VECSTANDARD */
2073: *a = *((PetscScalar **)x->data);
2074: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
2075: PetscFunctionReturn(PETSC_SUCCESS);
2076: }
2078: /*@C
2079: VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`
2081: Not Collective
2083: Input Parameters:
2084: + x - the vector
2085: - a - the array
2087: Level: beginner
2089: Fortran Notes:
2090: `VecRestoreArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayReadF90()`
2092: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2093: @*/
2094: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar **a)
2095: {
2096: PetscFunctionBegin;
2098: if (a) PetscAssertPointer(a, 2);
2099: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
2100: /* nothing */
2101: } else if (x->ops->restorearrayread) { /* VECNEST */
2102: PetscUseTypeMethod(x, restorearrayread, a);
2103: } else { /* No one? */
2104: PetscObjectState state;
2106: // x->ops->restorearray may bump the object state, but since we know this is a read-restore
2107: // we can just undo that
2108: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2109: PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
2110: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2111: }
2112: if (a) *a = NULL;
2113: PetscFunctionReturn(PETSC_SUCCESS);
2114: }
2116: /*@C
2117: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
2118: MPI processes's portion of the vector data.
2120: Logically Collective
2122: Input Parameter:
2123: . x - the vector
2125: Output Parameter:
2126: . a - location to put pointer to the array
2128: Level: intermediate
2130: Note:
2131: The values in this array are NOT valid, the caller of this routine is responsible for putting
2132: values into the array; any values it does not set will be invalid.
2134: The array must be returned using a matching call to `VecRestoreArrayWrite()`.
2136: For vectors associated with GPUs, the host and device vectors are not synchronized before
2137: giving access. If you need correct values in the array use `VecGetArray()`
2139: Fortran Notes:
2140: `VecGetArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayWriteF90()`
2142: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteF90()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
2143: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
2144: @*/
2145: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar **a)
2146: {
2147: PetscFunctionBegin;
2149: PetscAssertPointer(a, 2);
2150: PetscCall(VecSetErrorIfLocked(x, 1));
2151: if (x->ops->getarraywrite) {
2152: PetscUseTypeMethod(x, getarraywrite, a);
2153: } else {
2154: PetscCall(VecGetArray(x, a));
2155: }
2156: PetscFunctionReturn(PETSC_SUCCESS);
2157: }
2159: /*@C
2160: VecRestoreArrayWrite - Restores a vector after `VecGetArrayWrite()` has been called.
2162: Logically Collective
2164: Input Parameters:
2165: + x - the vector
2166: - a - location of pointer to array obtained from `VecGetArray()`
2168: Level: beginner
2170: Fortran Notes:
2171: `VecRestoreArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayWriteF90()`
2173: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteF90()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2174: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2175: @*/
2176: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar **a)
2177: {
2178: PetscFunctionBegin;
2180: if (a) PetscAssertPointer(a, 2);
2181: if (x->ops->restorearraywrite) {
2182: PetscUseTypeMethod(x, restorearraywrite, a);
2183: } else if (x->ops->restorearray) {
2184: PetscUseTypeMethod(x, restorearray, a);
2185: }
2186: if (a) *a = NULL;
2187: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2188: PetscFunctionReturn(PETSC_SUCCESS);
2189: }
2191: /*@C
2192: VecGetArrays - Returns a pointer to the arrays in a set of vectors
2193: that were created by a call to `VecDuplicateVecs()`.
2195: Logically Collective; No Fortran Support
2197: Input Parameters:
2198: + x - the vectors
2199: - n - the number of vectors
2201: Output Parameter:
2202: . a - location to put pointer to the array
2204: Level: intermediate
2206: Note:
2207: You MUST call `VecRestoreArrays()` when you no longer need access to the arrays.
2209: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2210: @*/
2211: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2212: {
2213: PetscInt i;
2214: PetscScalar **q;
2216: PetscFunctionBegin;
2217: PetscAssertPointer(x, 1);
2219: PetscAssertPointer(a, 3);
2220: PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2221: PetscCall(PetscMalloc1(n, &q));
2222: for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2223: *a = q;
2224: PetscFunctionReturn(PETSC_SUCCESS);
2225: }
2227: /*@C
2228: VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2229: has been called.
2231: Logically Collective; No Fortran Support
2233: Input Parameters:
2234: + x - the vector
2235: . n - the number of vectors
2236: - a - location of pointer to arrays obtained from `VecGetArrays()`
2238: Notes:
2239: For regular PETSc vectors this routine does not involve any copies. For
2240: any special vectors that do not store local vector data in a contiguous
2241: array, this routine will copy the data back into the underlying
2242: vector data structure from the arrays obtained with `VecGetArrays()`.
2244: Level: intermediate
2246: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2247: @*/
2248: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2249: {
2250: PetscInt i;
2251: PetscScalar **q = *a;
2253: PetscFunctionBegin;
2254: PetscAssertPointer(x, 1);
2256: PetscAssertPointer(a, 3);
2258: for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2259: PetscCall(PetscFree(q));
2260: PetscFunctionReturn(PETSC_SUCCESS);
2261: }
2263: /*@C
2264: VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g.,
2265: `VECCUDA`), the returned pointer will be a device pointer to the device memory that contains
2266: this MPI processes's portion of the vector data.
2268: Logically Collective; No Fortran Support
2270: Input Parameter:
2271: . x - the vector
2273: Output Parameters:
2274: + a - location to put pointer to the array
2275: - mtype - memory type of the array
2277: Level: beginner
2279: Note:
2280: Device data is guaranteed to have the latest value. Otherwise, when this is a host vector
2281: (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host
2282: pointer.
2284: For `VECKOKKOS`, if Kokkos is configured without device (e.g., use serial or openmp), per
2285: this function, the vector works like `VECSEQ`/`VECMPI`; otherwise, it works like `VECCUDA` or
2286: `VECHIP` etc.
2288: Use `VecRestoreArrayAndMemType()` when the array access is no longer needed.
2290: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`,
2291: `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2292: @*/
2293: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2294: {
2295: PetscFunctionBegin;
2298: PetscAssertPointer(a, 2);
2299: if (mtype) PetscAssertPointer(mtype, 3);
2300: PetscCall(VecSetErrorIfLocked(x, 1));
2301: if (x->ops->getarrayandmemtype) {
2302: /* VECCUDA, VECKOKKOS etc */
2303: PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2304: } else {
2305: /* VECSTANDARD, VECNEST, VECVIENNACL */
2306: PetscCall(VecGetArray(x, a));
2307: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2308: }
2309: PetscFunctionReturn(PETSC_SUCCESS);
2310: }
2312: /*@C
2313: VecRestoreArrayAndMemType - Restores a vector after `VecGetArrayAndMemType()` has been called.
2315: Logically Collective; No Fortran Support
2317: Input Parameters:
2318: + x - the vector
2319: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`
2321: Level: beginner
2323: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`,
2324: `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2325: @*/
2326: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar **a)
2327: {
2328: PetscFunctionBegin;
2331: if (a) PetscAssertPointer(a, 2);
2332: if (x->ops->restorearrayandmemtype) {
2333: /* VECCUDA, VECKOKKOS etc */
2334: PetscUseTypeMethod(x, restorearrayandmemtype, a);
2335: } else {
2336: /* VECNEST, VECVIENNACL */
2337: PetscCall(VecRestoreArray(x, a));
2338: } /* VECSTANDARD does nothing */
2339: if (a) *a = NULL;
2340: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2341: PetscFunctionReturn(PETSC_SUCCESS);
2342: }
2344: /*@C
2345: VecGetArrayReadAndMemType - Like `VecGetArrayRead()`, but if the input vector is a device vector, it will return a read-only device pointer.
2346: The returned pointer is guaranteed to point to up-to-date data. For host vectors, it functions as `VecGetArrayRead()`.
2348: Not Collective; No Fortran Support
2350: Input Parameter:
2351: . x - the vector
2353: Output Parameters:
2354: + a - the array
2355: - mtype - memory type of the array
2357: Level: beginner
2359: Notes:
2360: The array must be returned using a matching call to `VecRestoreArrayReadAndMemType()`.
2362: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2363: @*/
2364: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar **a, PetscMemType *mtype)
2365: {
2366: PetscFunctionBegin;
2369: PetscAssertPointer(a, 2);
2370: if (mtype) PetscAssertPointer(mtype, 3);
2371: if (x->ops->getarrayreadandmemtype) {
2372: /* VECCUDA/VECHIP though they are also petscnative */
2373: PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2374: } else if (x->ops->getarrayandmemtype) {
2375: /* VECKOKKOS */
2376: PetscObjectState state;
2378: // see VecGetArrayRead() for why
2379: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2380: PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2381: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2382: } else {
2383: PetscCall(VecGetArrayRead(x, a));
2384: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2385: }
2386: PetscFunctionReturn(PETSC_SUCCESS);
2387: }
2389: /*@C
2390: VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`
2392: Not Collective; No Fortran Support
2394: Input Parameters:
2395: + x - the vector
2396: - a - the array
2398: Level: beginner
2400: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2401: @*/
2402: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar **a)
2403: {
2404: PetscFunctionBegin;
2407: if (a) PetscAssertPointer(a, 2);
2408: if (x->ops->restorearrayreadandmemtype) {
2409: /* VECCUDA/VECHIP */
2410: PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2411: } else if (!x->petscnative) {
2412: /* VECNEST */
2413: PetscCall(VecRestoreArrayRead(x, a));
2414: }
2415: if (a) *a = NULL;
2416: PetscFunctionReturn(PETSC_SUCCESS);
2417: }
2419: /*@C
2420: VecGetArrayWriteAndMemType - Like `VecGetArrayWrite()`, but if this is a device vector it will always return
2421: a device pointer to the device memory that contains this processor's portion of the vector data.
2423: Logically Collective; No Fortran Support
2425: Input Parameter:
2426: . x - the vector
2428: Output Parameters:
2429: + a - the array
2430: - mtype - memory type of the array
2432: Level: beginner
2434: Note:
2435: The array must be returned using a matching call to `VecRestoreArrayWriteAndMemType()`, where it will label the device memory as most recent.
2437: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2438: @*/
2439: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2440: {
2441: PetscFunctionBegin;
2444: PetscCall(VecSetErrorIfLocked(x, 1));
2445: PetscAssertPointer(a, 2);
2446: if (mtype) PetscAssertPointer(mtype, 3);
2447: if (x->ops->getarraywriteandmemtype) {
2448: /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2449: PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2450: } else if (x->ops->getarrayandmemtype) {
2451: PetscCall(VecGetArrayAndMemType(x, a, mtype));
2452: } else {
2453: /* VECNEST, VECVIENNACL */
2454: PetscCall(VecGetArrayWrite(x, a));
2455: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2456: }
2457: PetscFunctionReturn(PETSC_SUCCESS);
2458: }
2460: /*@C
2461: VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`
2463: Logically Collective; No Fortran Support
2465: Input Parameters:
2466: + x - the vector
2467: - a - the array
2469: Level: beginner
2471: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2472: @*/
2473: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar **a)
2474: {
2475: PetscFunctionBegin;
2478: PetscCall(VecSetErrorIfLocked(x, 1));
2479: if (a) PetscAssertPointer(a, 2);
2480: if (x->ops->restorearraywriteandmemtype) {
2481: /* VECCUDA/VECHIP */
2482: PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2483: PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2484: } else if (x->ops->restorearrayandmemtype) {
2485: PetscCall(VecRestoreArrayAndMemType(x, a));
2486: } else {
2487: PetscCall(VecRestoreArray(x, a));
2488: }
2489: if (a) *a = NULL;
2490: PetscFunctionReturn(PETSC_SUCCESS);
2491: }
2493: /*@
2494: VecPlaceArray - Allows one to replace the array in a vector with an
2495: array provided by the user. This is useful to avoid copying an array
2496: into a vector.
2498: Logically Collective; No Fortran Support
2500: Input Parameters:
2501: + vec - the vector
2502: - array - the array
2504: Level: developer
2506: Notes:
2507: Use `VecReplaceArray()` instead to permanently replace the array
2509: You can return to the original array with a call to `VecResetArray()`. `vec` does not take
2510: ownership of `array` in any way.
2512: The user must free `array` themselves but be careful not to
2513: do so before the vector has either been destroyed, had its original array restored with
2514: `VecResetArray()` or permanently replaced with `VecReplaceArray()`.
2516: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2517: @*/
2518: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2519: {
2520: PetscFunctionBegin;
2523: if (array) PetscAssertPointer(array, 2);
2524: PetscUseTypeMethod(vec, placearray, array);
2525: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2526: PetscFunctionReturn(PETSC_SUCCESS);
2527: }
2529: /*@C
2530: VecReplaceArray - 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; No Fortran Support
2536: Input Parameters:
2537: + vec - the vector
2538: - array - the array
2540: Level: developer
2542: Notes:
2543: This permanently replaces the array and frees the memory associated
2544: with the old array. Use `VecPlaceArray()` to temporarily replace the array.
2546: The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
2547: freed by the user. It will be freed when the vector is destroyed.
2549: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2550: @*/
2551: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2552: {
2553: PetscFunctionBegin;
2556: PetscUseTypeMethod(vec, replacearray, array);
2557: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2558: PetscFunctionReturn(PETSC_SUCCESS);
2559: }
2561: /*MC
2562: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
2563: and makes them accessible via a Fortran pointer.
2565: Synopsis:
2566: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
2568: Collective
2570: Input Parameters:
2571: + x - a vector to mimic
2572: - n - the number of vectors to obtain
2574: Output Parameters:
2575: + y - Fortran pointer to the array of vectors
2576: - ierr - error code
2578: Example of Usage:
2579: .vb
2580: #include <petsc/finclude/petscvec.h>
2581: use petscvec
2583: Vec x
2584: Vec, pointer :: y(:)
2585: ....
2586: call VecDuplicateVecsF90(x,2,y,ierr)
2587: call VecSet(y(2),alpha,ierr)
2588: call VecSet(y(2),alpha,ierr)
2589: ....
2590: call VecDestroyVecsF90(2,y,ierr)
2591: .ve
2593: Level: beginner
2595: Note:
2596: Use `VecDestroyVecsF90()` to free the space.
2598: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecsF90()`, `VecDuplicateVecs()`
2599: M*/
2601: /*MC
2602: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2603: `VecGetArrayF90()`.
2605: Synopsis:
2606: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2608: Logically Collective
2610: Input Parameters:
2611: + x - vector
2612: - xx_v - the Fortran pointer to the array
2614: Output Parameter:
2615: . ierr - error code
2617: Example of Usage:
2618: .vb
2619: #include <petsc/finclude/petscvec.h>
2620: use petscvec
2622: PetscScalar, pointer :: xx_v(:)
2623: ....
2624: call VecGetArrayF90(x,xx_v,ierr)
2625: xx_v(3) = a
2626: call VecRestoreArrayF90(x,xx_v,ierr)
2627: .ve
2629: Level: beginner
2631: .seealso: [](ch_vectors), `Vec`, `VecGetArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrayReadF90()`
2632: M*/
2634: /*MC
2635: VecDestroyVecsF90 - Frees a block of vectors obtained with `VecDuplicateVecsF90()`.
2637: Synopsis:
2638: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
2640: Collective
2642: Input Parameters:
2643: + n - the number of vectors previously obtained
2644: - x - pointer to array of vector pointers
2646: Output Parameter:
2647: . ierr - error code
2649: Level: beginner
2651: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecs()`, `VecDuplicateVecsF90()`
2652: M*/
2654: /*MC
2655: VecGetArrayF90 - Accesses a vector array from Fortran. For default PETSc
2656: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2657: this routine is implementation dependent. You MUST call `VecRestoreArrayF90()`
2658: when you no longer need access to the array.
2660: Synopsis:
2661: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2663: Logically Collective
2665: Input Parameter:
2666: . x - vector
2668: Output Parameters:
2669: + xx_v - the Fortran pointer to the array
2670: - ierr - error code
2672: Example of Usage:
2673: .vb
2674: #include <petsc/finclude/petscvec.h>
2675: use petscvec
2677: PetscScalar, pointer :: xx_v(:)
2678: ....
2679: call VecGetArrayF90(x,xx_v,ierr)
2680: xx_v(3) = a
2681: call VecRestoreArrayF90(x,xx_v,ierr)
2682: .ve
2684: Level: beginner
2686: Note:
2687: If you ONLY intend to read entries from the array and not change any entries you should use `VecGetArrayReadF90()`.
2689: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayReadF90()`
2690: M*/
2692: /*MC
2693: VecGetArrayReadF90 - Accesses a read only array from Fortran. For default PETSc
2694: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2695: this routine is implementation dependent. You MUST call `VecRestoreArrayReadF90()`
2696: when you no longer need access to the array.
2698: Synopsis:
2699: VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2701: Logically Collective
2703: Input Parameter:
2704: . x - vector
2706: Output Parameters:
2707: + xx_v - the Fortran pointer to the array
2708: - ierr - error code
2710: Example of Usage:
2711: .vb
2712: #include <petsc/finclude/petscvec.h>
2713: use petscvec
2715: PetscScalar, pointer :: xx_v(:)
2716: ....
2717: call VecGetArrayReadF90(x,xx_v,ierr)
2718: a = xx_v(3)
2719: call VecRestoreArrayReadF90(x,xx_v,ierr)
2720: .ve
2722: Level: beginner
2724: Note:
2725: If you intend to write entries into the array you must use `VecGetArrayF90()`.
2727: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecGetArrayF90()`
2728: M*/
2730: /*MC
2731: VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2732: `VecGetArrayReadF90()`.
2734: Synopsis:
2735: VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2737: Logically Collective
2739: Input Parameters:
2740: + x - vector
2741: - xx_v - the Fortran pointer to the array
2743: Output Parameter:
2744: . ierr - error code
2746: Example of Usage:
2747: .vb
2748: #include <petsc/finclude/petscvec.h>
2749: use petscvec
2751: PetscScalar, pointer :: xx_v(:)
2752: ....
2753: call VecGetArrayReadF90(x,xx_v,ierr)
2754: a = xx_v(3)
2755: call VecRestoreArrayReadF90(x,xx_v,ierr)
2756: .ve
2758: Level: beginner
2760: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecRestoreArrayF90()`
2761: M*/
2763: /*@C
2764: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2765: processor's portion of the vector data. You MUST call `VecRestoreArray2d()`
2766: when you no longer need access to the array.
2768: Logically Collective
2770: Input Parameters:
2771: + x - the vector
2772: . m - first dimension of two dimensional array
2773: . n - second dimension of two dimensional array
2774: . mstart - first index you will use in first coordinate direction (often 0)
2775: - nstart - first index in the second coordinate direction (often 0)
2777: Output Parameter:
2778: . a - location to put pointer to the array
2780: Level: developer
2782: Notes:
2783: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2784: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2785: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2786: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2788: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2790: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2791: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2792: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2793: @*/
2794: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2795: {
2796: PetscInt i, N;
2797: PetscScalar *aa;
2799: PetscFunctionBegin;
2801: PetscAssertPointer(a, 6);
2803: PetscCall(VecGetLocalSize(x, &N));
2804: 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);
2805: PetscCall(VecGetArray(x, &aa));
2807: PetscCall(PetscMalloc1(m, a));
2808: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2809: *a -= mstart;
2810: PetscFunctionReturn(PETSC_SUCCESS);
2811: }
2813: /*@C
2814: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2815: processor's portion of the vector data. You MUST call `VecRestoreArray2dWrite()`
2816: when you no longer need access to the array.
2818: Logically Collective
2820: Input Parameters:
2821: + x - the vector
2822: . m - first dimension of two dimensional array
2823: . n - second dimension of two dimensional array
2824: . mstart - first index you will use in first coordinate direction (often 0)
2825: - nstart - first index in the second coordinate direction (often 0)
2827: Output Parameter:
2828: . a - location to put pointer to the array
2830: Level: developer
2832: Notes:
2833: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2834: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2835: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2836: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2838: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2840: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2841: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2842: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2843: @*/
2844: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2845: {
2846: PetscInt i, N;
2847: PetscScalar *aa;
2849: PetscFunctionBegin;
2851: PetscAssertPointer(a, 6);
2853: PetscCall(VecGetLocalSize(x, &N));
2854: 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);
2855: PetscCall(VecGetArrayWrite(x, &aa));
2857: PetscCall(PetscMalloc1(m, a));
2858: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2859: *a -= mstart;
2860: PetscFunctionReturn(PETSC_SUCCESS);
2861: }
2863: /*@C
2864: VecRestoreArray2d - Restores a vector after `VecGetArray2d()` has been called.
2866: Logically Collective
2868: Input Parameters:
2869: + x - the vector
2870: . m - first dimension of two dimensional array
2871: . n - second dimension of the two dimensional array
2872: . mstart - first index you will use in first coordinate direction (often 0)
2873: . nstart - first index in the second coordinate direction (often 0)
2874: - a - location of pointer to array obtained from `VecGetArray2d()`
2876: Level: developer
2878: Notes:
2879: For regular PETSc vectors this routine does not involve any copies. For
2880: any special vectors that do not store local vector data in a contiguous
2881: array, this routine will copy the data back into the underlying
2882: vector data structure from the array obtained with `VecGetArray()`.
2884: This routine actually zeros out the `a` pointer.
2886: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2887: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2888: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2889: @*/
2890: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2891: {
2892: void *dummy;
2894: PetscFunctionBegin;
2896: PetscAssertPointer(a, 6);
2898: dummy = (void *)(*a + mstart);
2899: PetscCall(PetscFree(dummy));
2900: PetscCall(VecRestoreArray(x, NULL));
2901: *a = NULL;
2902: PetscFunctionReturn(PETSC_SUCCESS);
2903: }
2905: /*@C
2906: VecRestoreArray2dWrite - Restores a vector after `VecGetArray2dWrite()` has been called.
2908: Logically Collective
2910: Input Parameters:
2911: + x - the vector
2912: . m - first dimension of two dimensional array
2913: . n - second dimension of the two dimensional array
2914: . mstart - first index you will use in first coordinate direction (often 0)
2915: . nstart - first index in the second coordinate direction (often 0)
2916: - a - location of pointer to array obtained from `VecGetArray2d()`
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 `VecGetArray()`.
2926: This routine actually zeros out the `a` pointer.
2928: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2929: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2930: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2931: @*/
2932: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2933: {
2934: void *dummy;
2936: PetscFunctionBegin;
2938: PetscAssertPointer(a, 6);
2940: dummy = (void *)(*a + mstart);
2941: PetscCall(PetscFree(dummy));
2942: PetscCall(VecRestoreArrayWrite(x, NULL));
2943: PetscFunctionReturn(PETSC_SUCCESS);
2944: }
2946: /*@C
2947: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2948: processor's portion of the vector data. You MUST call `VecRestoreArray1d()`
2949: when you no longer need access to the array.
2951: Logically Collective
2953: Input Parameters:
2954: + x - the vector
2955: . m - first dimension of two dimensional array
2956: - mstart - first index you will use in first 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` is likely
2965: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2966: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2968: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2970: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2971: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2972: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2973: @*/
2974: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2975: {
2976: PetscInt N;
2978: PetscFunctionBegin;
2980: PetscAssertPointer(a, 4);
2982: PetscCall(VecGetLocalSize(x, &N));
2983: 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);
2984: PetscCall(VecGetArray(x, a));
2985: *a -= mstart;
2986: PetscFunctionReturn(PETSC_SUCCESS);
2987: }
2989: /*@C
2990: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
2991: processor's portion of the vector data. You MUST call `VecRestoreArray1dWrite()`
2992: when you no longer need access to the array.
2994: Logically Collective
2996: Input Parameters:
2997: + x - the vector
2998: . m - first dimension of two dimensional array
2999: - mstart - first index you will use in first coordinate direction (often 0)
3001: Output Parameter:
3002: . a - location to put pointer to the array
3004: Level: developer
3006: Notes:
3007: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3008: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3009: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3011: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3013: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3014: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3015: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3016: @*/
3017: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3018: {
3019: PetscInt N;
3021: PetscFunctionBegin;
3023: PetscAssertPointer(a, 4);
3025: PetscCall(VecGetLocalSize(x, &N));
3026: 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);
3027: PetscCall(VecGetArrayWrite(x, a));
3028: *a -= mstart;
3029: PetscFunctionReturn(PETSC_SUCCESS);
3030: }
3032: /*@C
3033: VecRestoreArray1d - Restores a vector after `VecGetArray1d()` has been called.
3035: Logically Collective
3037: Input Parameters:
3038: + x - the vector
3039: . m - first dimension of two dimensional array
3040: . mstart - first index you will use in first coordinate direction (often 0)
3041: - a - location of pointer to array obtained from `VecGetArray1d()`
3043: Level: developer
3045: Notes:
3046: For regular PETSc vectors this routine does not involve any copies. For
3047: any special vectors that do not store local vector data in a contiguous
3048: array, this routine will copy the data back into the underlying
3049: vector data structure from the array obtained with `VecGetArray1d()`.
3051: This routine actually zeros out the `a` pointer.
3053: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3054: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3055: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3056: @*/
3057: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3058: {
3059: PetscFunctionBegin;
3062: PetscCall(VecRestoreArray(x, NULL));
3063: *a = NULL;
3064: PetscFunctionReturn(PETSC_SUCCESS);
3065: }
3067: /*@C
3068: VecRestoreArray1dWrite - Restores a vector after `VecGetArray1dWrite()` has been called.
3070: Logically Collective
3072: Input Parameters:
3073: + x - the vector
3074: . m - first dimension of two dimensional array
3075: . mstart - first index you will use in first coordinate direction (often 0)
3076: - a - location of pointer to array obtained from `VecGetArray1d()`
3078: Level: developer
3080: Notes:
3081: For regular PETSc vectors this routine does not involve any copies. For
3082: any special vectors that do not store local vector data in a contiguous
3083: array, this routine will copy the data back into the underlying
3084: vector data structure from the array obtained with `VecGetArray1d()`.
3086: This routine actually zeros out the `a` pointer.
3088: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3089: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3090: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3091: @*/
3092: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3093: {
3094: PetscFunctionBegin;
3097: PetscCall(VecRestoreArrayWrite(x, NULL));
3098: *a = NULL;
3099: PetscFunctionReturn(PETSC_SUCCESS);
3100: }
3102: /*@C
3103: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
3104: processor's portion of the vector data. You MUST call `VecRestoreArray3d()`
3105: when you no longer need access to the array.
3107: Logically Collective
3109: Input Parameters:
3110: + x - the vector
3111: . m - first dimension of three dimensional array
3112: . n - second dimension of three dimensional array
3113: . p - third dimension of three dimensional array
3114: . mstart - first index you will use in first coordinate direction (often 0)
3115: . nstart - first index in the second coordinate direction (often 0)
3116: - pstart - first index in the third coordinate direction (often 0)
3118: Output Parameter:
3119: . a - location to put pointer to the array
3121: Level: developer
3123: Notes:
3124: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3125: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3126: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3127: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3129: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3131: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3132: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecRestoreArray3d()`,
3133: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3134: @*/
3135: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3136: {
3137: PetscInt i, N, j;
3138: PetscScalar *aa, **b;
3140: PetscFunctionBegin;
3142: PetscAssertPointer(a, 8);
3144: PetscCall(VecGetLocalSize(x, &N));
3145: 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);
3146: PetscCall(VecGetArray(x, &aa));
3148: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3149: b = (PetscScalar **)((*a) + m);
3150: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3151: for (i = 0; i < m; i++)
3152: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset(aa, i * n * p + j * p - pstart);
3153: *a -= mstart;
3154: PetscFunctionReturn(PETSC_SUCCESS);
3155: }
3157: /*@C
3158: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
3159: processor's portion of the vector data. You MUST call `VecRestoreArray3dWrite()`
3160: when you no longer need access to the array.
3162: Logically Collective
3164: Input Parameters:
3165: + x - the vector
3166: . m - first dimension of three dimensional array
3167: . n - second dimension of three dimensional array
3168: . p - third dimension of three dimensional array
3169: . mstart - first index you will use in first coordinate direction (often 0)
3170: . nstart - first index in the second coordinate direction (often 0)
3171: - pstart - first index in the third coordinate direction (often 0)
3173: Output Parameter:
3174: . a - location to put pointer to the array
3176: Level: developer
3178: Notes:
3179: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3180: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3181: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3182: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3184: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3186: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3187: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3188: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3189: @*/
3190: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3191: {
3192: PetscInt i, N, j;
3193: PetscScalar *aa, **b;
3195: PetscFunctionBegin;
3197: PetscAssertPointer(a, 8);
3199: PetscCall(VecGetLocalSize(x, &N));
3200: 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);
3201: PetscCall(VecGetArrayWrite(x, &aa));
3203: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3204: b = (PetscScalar **)((*a) + m);
3205: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3206: for (i = 0; i < m; i++)
3207: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3209: *a -= mstart;
3210: PetscFunctionReturn(PETSC_SUCCESS);
3211: }
3213: /*@C
3214: VecRestoreArray3d - Restores a vector after `VecGetArray3d()` has been called.
3216: Logically Collective
3218: Input Parameters:
3219: + x - the vector
3220: . m - first dimension of three dimensional array
3221: . n - second dimension of the three dimensional array
3222: . p - third dimension of the three dimensional array
3223: . mstart - first index you will use in first coordinate direction (often 0)
3224: . nstart - first index in the second coordinate direction (often 0)
3225: . pstart - first index in the third coordinate direction (often 0)
3226: - a - location of pointer to array obtained from VecGetArray3d()
3228: Level: developer
3230: Notes:
3231: For regular PETSc vectors this routine does not involve any copies. For
3232: any special vectors that do not store local vector data in a contiguous
3233: array, this routine will copy the data back into the underlying
3234: vector data structure from the array obtained with `VecGetArray()`.
3236: This routine actually zeros out the `a` pointer.
3238: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3239: `VecGetArray2d()`, `VecGetArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3240: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3241: @*/
3242: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3243: {
3244: void *dummy;
3246: PetscFunctionBegin;
3248: PetscAssertPointer(a, 8);
3250: dummy = (void *)(*a + mstart);
3251: PetscCall(PetscFree(dummy));
3252: PetscCall(VecRestoreArray(x, NULL));
3253: *a = NULL;
3254: PetscFunctionReturn(PETSC_SUCCESS);
3255: }
3257: /*@C
3258: VecRestoreArray3dWrite - Restores a vector after `VecGetArray3dWrite()` has been called.
3260: Logically Collective
3262: Input Parameters:
3263: + x - the vector
3264: . m - first dimension of three dimensional array
3265: . n - second dimension of the three dimensional array
3266: . p - third dimension of the three dimensional array
3267: . mstart - first index you will use in first coordinate direction (often 0)
3268: . nstart - first index in the second coordinate direction (often 0)
3269: . pstart - first index in the third coordinate direction (often 0)
3270: - a - location of pointer to array obtained from VecGetArray3d()
3272: Level: developer
3274: Notes:
3275: For regular PETSc vectors this routine does not involve any copies. For
3276: any special vectors that do not store local vector data in a contiguous
3277: array, this routine will copy the data back into the underlying
3278: vector data structure from the array obtained with `VecGetArray()`.
3280: This routine actually zeros out the `a` pointer.
3282: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3283: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3284: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3285: @*/
3286: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3287: {
3288: void *dummy;
3290: PetscFunctionBegin;
3292: PetscAssertPointer(a, 8);
3294: dummy = (void *)(*a + mstart);
3295: PetscCall(PetscFree(dummy));
3296: PetscCall(VecRestoreArrayWrite(x, NULL));
3297: *a = NULL;
3298: PetscFunctionReturn(PETSC_SUCCESS);
3299: }
3301: /*@C
3302: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3303: processor's portion of the vector data. You MUST call `VecRestoreArray4d()`
3304: when you no longer need access to the array.
3306: Logically Collective
3308: Input Parameters:
3309: + x - the vector
3310: . m - first dimension of four dimensional array
3311: . n - second dimension of four dimensional array
3312: . p - third dimension of four dimensional array
3313: . q - fourth dimension of four dimensional array
3314: . mstart - first index you will use in first coordinate direction (often 0)
3315: . nstart - first index in the second coordinate direction (often 0)
3316: . pstart - first index in the third coordinate direction (often 0)
3317: - qstart - first index in the fourth coordinate direction (often 0)
3319: Output Parameter:
3320: . a - location to put pointer to the array
3322: Level: developer
3324: Notes:
3325: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3326: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3327: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3328: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3330: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3332: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3333: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3334: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecRestoreArray4d()`
3335: @*/
3336: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3337: {
3338: PetscInt i, N, j, k;
3339: PetscScalar *aa, ***b, **c;
3341: PetscFunctionBegin;
3343: PetscAssertPointer(a, 10);
3345: PetscCall(VecGetLocalSize(x, &N));
3346: 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);
3347: PetscCall(VecGetArray(x, &aa));
3349: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3350: b = (PetscScalar ***)((*a) + m);
3351: c = (PetscScalar **)(b + m * n);
3352: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3353: for (i = 0; i < m; i++)
3354: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3355: for (i = 0; i < m; i++)
3356: for (j = 0; j < n; j++)
3357: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3358: *a -= mstart;
3359: PetscFunctionReturn(PETSC_SUCCESS);
3360: }
3362: /*@C
3363: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3364: processor's portion of the vector data. You MUST call `VecRestoreArray4dWrite()`
3365: when you no longer need access to the array.
3367: Logically Collective
3369: Input Parameters:
3370: + x - the vector
3371: . m - first dimension of four dimensional array
3372: . n - second dimension of four dimensional array
3373: . p - third dimension of four dimensional array
3374: . q - fourth dimension of four dimensional array
3375: . mstart - first index you will use in first coordinate direction (often 0)
3376: . nstart - first index in the second coordinate direction (often 0)
3377: . pstart - first index in the third coordinate direction (often 0)
3378: - qstart - first index in the fourth coordinate direction (often 0)
3380: Output Parameter:
3381: . a - location to put pointer to the array
3383: Level: developer
3385: Notes:
3386: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3387: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3388: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3389: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3391: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3393: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3394: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3395: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3396: @*/
3397: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3398: {
3399: PetscInt i, N, j, k;
3400: PetscScalar *aa, ***b, **c;
3402: PetscFunctionBegin;
3404: PetscAssertPointer(a, 10);
3406: PetscCall(VecGetLocalSize(x, &N));
3407: 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);
3408: PetscCall(VecGetArrayWrite(x, &aa));
3410: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3411: b = (PetscScalar ***)((*a) + m);
3412: c = (PetscScalar **)(b + m * n);
3413: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3414: for (i = 0; i < m; i++)
3415: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3416: for (i = 0; i < m; i++)
3417: for (j = 0; j < n; j++)
3418: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3419: *a -= mstart;
3420: PetscFunctionReturn(PETSC_SUCCESS);
3421: }
3423: /*@C
3424: VecRestoreArray4d - Restores a vector after `VecGetArray4d()` has been called.
3426: Logically Collective
3428: Input Parameters:
3429: + x - the vector
3430: . m - first dimension of four dimensional array
3431: . n - second dimension of the four dimensional array
3432: . p - third dimension of the four dimensional array
3433: . q - fourth dimension of the four dimensional array
3434: . mstart - first index you will use in first coordinate direction (often 0)
3435: . nstart - first index in the second coordinate direction (often 0)
3436: . pstart - first index in the third coordinate direction (often 0)
3437: . qstart - first index in the fourth coordinate direction (often 0)
3438: - a - location of pointer to array obtained from VecGetArray4d()
3440: Level: developer
3442: Notes:
3443: For regular PETSc vectors this routine does not involve any copies. For
3444: any special vectors that do not store local vector data in a contiguous
3445: array, this routine will copy the data back into the underlying
3446: vector data structure from the array obtained with `VecGetArray()`.
3448: This routine actually zeros out the `a` pointer.
3450: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3451: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3452: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecGet`
3453: @*/
3454: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3455: {
3456: void *dummy;
3458: PetscFunctionBegin;
3460: PetscAssertPointer(a, 10);
3462: dummy = (void *)(*a + mstart);
3463: PetscCall(PetscFree(dummy));
3464: PetscCall(VecRestoreArray(x, NULL));
3465: *a = NULL;
3466: PetscFunctionReturn(PETSC_SUCCESS);
3467: }
3469: /*@C
3470: VecRestoreArray4dWrite - Restores a vector after `VecGetArray4dWrite()` has been called.
3472: Logically Collective
3474: Input Parameters:
3475: + x - the vector
3476: . m - first dimension of four dimensional array
3477: . n - second dimension of the four dimensional array
3478: . p - third dimension of the four dimensional array
3479: . q - fourth dimension of the four dimensional array
3480: . mstart - first index you will use in first coordinate direction (often 0)
3481: . nstart - first index in the second coordinate direction (often 0)
3482: . pstart - first index in the third coordinate direction (often 0)
3483: . qstart - first index in the fourth coordinate direction (often 0)
3484: - a - location of pointer to array obtained from `VecGetArray4d()`
3486: Level: developer
3488: Notes:
3489: For regular PETSc vectors this routine does not involve any copies. For
3490: any special vectors that do not store local vector data in a contiguous
3491: array, this routine will copy the data back into the underlying
3492: vector data structure from the array obtained with `VecGetArray()`.
3494: This routine actually zeros out the `a` pointer.
3496: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3497: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3498: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3499: @*/
3500: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3501: {
3502: void *dummy;
3504: PetscFunctionBegin;
3506: PetscAssertPointer(a, 10);
3508: dummy = (void *)(*a + mstart);
3509: PetscCall(PetscFree(dummy));
3510: PetscCall(VecRestoreArrayWrite(x, NULL));
3511: *a = NULL;
3512: PetscFunctionReturn(PETSC_SUCCESS);
3513: }
3515: /*@C
3516: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3517: processor's portion of the vector data. You MUST call `VecRestoreArray2dRead()`
3518: when you no longer need access to the array.
3520: Logically Collective
3522: Input Parameters:
3523: + x - the vector
3524: . m - first dimension of two dimensional array
3525: . n - second dimension of two dimensional array
3526: . mstart - first index you will use in first coordinate direction (often 0)
3527: - nstart - first index in the second coordinate direction (often 0)
3529: Output Parameter:
3530: . a - location to put pointer to the array
3532: Level: developer
3534: Notes:
3535: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
3536: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3537: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3538: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
3540: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3542: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3543: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3544: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3545: @*/
3546: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3547: {
3548: PetscInt i, N;
3549: const PetscScalar *aa;
3551: PetscFunctionBegin;
3553: PetscAssertPointer(a, 6);
3555: PetscCall(VecGetLocalSize(x, &N));
3556: 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);
3557: PetscCall(VecGetArrayRead(x, &aa));
3559: PetscCall(PetscMalloc1(m, a));
3560: for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3561: *a -= mstart;
3562: PetscFunctionReturn(PETSC_SUCCESS);
3563: }
3565: /*@C
3566: VecRestoreArray2dRead - Restores a vector after `VecGetArray2dRead()` has been called.
3568: Logically Collective
3570: Input Parameters:
3571: + x - the vector
3572: . m - first dimension of two dimensional array
3573: . n - second dimension of the two dimensional array
3574: . mstart - first index you will use in first coordinate direction (often 0)
3575: . nstart - first index in the second coordinate direction (often 0)
3576: - a - location of pointer to array obtained from VecGetArray2d()
3578: Level: developer
3580: Notes:
3581: For regular PETSc vectors this routine does not involve any copies. For
3582: any special vectors that do not store local vector data in a contiguous
3583: array, this routine will copy the data back into the underlying
3584: vector data structure from the array obtained with `VecGetArray()`.
3586: This routine actually zeros out the `a` pointer.
3588: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3589: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3590: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3591: @*/
3592: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3593: {
3594: void *dummy;
3596: PetscFunctionBegin;
3598: PetscAssertPointer(a, 6);
3600: dummy = (void *)(*a + mstart);
3601: PetscCall(PetscFree(dummy));
3602: PetscCall(VecRestoreArrayRead(x, NULL));
3603: *a = NULL;
3604: PetscFunctionReturn(PETSC_SUCCESS);
3605: }
3607: /*@C
3608: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3609: processor's portion of the vector data. You MUST call `VecRestoreArray1dRead()`
3610: when you no longer need access to the array.
3612: Logically Collective
3614: Input Parameters:
3615: + x - the vector
3616: . m - first dimension of two dimensional array
3617: - mstart - first index you will use in first coordinate direction (often 0)
3619: Output Parameter:
3620: . a - location to put pointer to the array
3622: Level: developer
3624: Notes:
3625: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3626: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3627: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3629: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3631: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3632: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3633: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3634: @*/
3635: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3636: {
3637: PetscInt N;
3639: PetscFunctionBegin;
3641: PetscAssertPointer(a, 4);
3643: PetscCall(VecGetLocalSize(x, &N));
3644: 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);
3645: PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3646: *a -= mstart;
3647: PetscFunctionReturn(PETSC_SUCCESS);
3648: }
3650: /*@C
3651: VecRestoreArray1dRead - Restores a vector after `VecGetArray1dRead()` has been called.
3653: Logically Collective
3655: Input Parameters:
3656: + x - the vector
3657: . m - first dimension of two dimensional array
3658: . mstart - first index you will use in first coordinate direction (often 0)
3659: - a - location of pointer to array obtained from `VecGetArray1dRead()`
3661: Level: developer
3663: Notes:
3664: For regular PETSc vectors this routine does not involve any copies. For
3665: any special vectors that do not store local vector data in a contiguous
3666: array, this routine will copy the data back into the underlying
3667: vector data structure from the array obtained with `VecGetArray1dRead()`.
3669: This routine actually zeros out the `a` pointer.
3671: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3672: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3673: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3674: @*/
3675: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3676: {
3677: PetscFunctionBegin;
3680: PetscCall(VecRestoreArrayRead(x, NULL));
3681: *a = NULL;
3682: PetscFunctionReturn(PETSC_SUCCESS);
3683: }
3685: /*@C
3686: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3687: processor's portion of the vector data. You MUST call `VecRestoreArray3dRead()`
3688: when you no longer need access to the array.
3690: Logically Collective
3692: Input Parameters:
3693: + x - the vector
3694: . m - first dimension of three dimensional array
3695: . n - second dimension of three dimensional array
3696: . p - third dimension of three dimensional array
3697: . mstart - first index you will use in first coordinate direction (often 0)
3698: . nstart - first index in the second coordinate direction (often 0)
3699: - pstart - first index in the third coordinate direction (often 0)
3701: Output Parameter:
3702: . a - location to put pointer to the array
3704: Level: developer
3706: Notes:
3707: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3708: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3709: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3710: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3dRead()`.
3712: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3714: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3715: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3716: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3717: @*/
3718: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3719: {
3720: PetscInt i, N, j;
3721: const PetscScalar *aa;
3722: PetscScalar **b;
3724: PetscFunctionBegin;
3726: PetscAssertPointer(a, 8);
3728: PetscCall(VecGetLocalSize(x, &N));
3729: 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);
3730: PetscCall(VecGetArrayRead(x, &aa));
3732: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3733: b = (PetscScalar **)((*a) + m);
3734: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3735: for (i = 0; i < m; i++)
3736: for (j = 0; j < n; j++) b[i * n + j] = PetscSafePointerPlusOffset((PetscScalar *)aa, i * n * p + j * p - pstart);
3737: *a -= mstart;
3738: PetscFunctionReturn(PETSC_SUCCESS);
3739: }
3741: /*@C
3742: VecRestoreArray3dRead - Restores a vector after `VecGetArray3dRead()` has been called.
3744: Logically Collective
3746: Input Parameters:
3747: + x - the vector
3748: . m - first dimension of three dimensional array
3749: . n - second dimension of the three dimensional array
3750: . p - third dimension of the three dimensional array
3751: . mstart - first index you will use in first coordinate direction (often 0)
3752: . nstart - first index in the second coordinate direction (often 0)
3753: . pstart - first index in the third coordinate direction (often 0)
3754: - a - location of pointer to array obtained from `VecGetArray3dRead()`
3756: Level: developer
3758: Notes:
3759: For regular PETSc vectors this routine does not involve any copies. For
3760: any special vectors that do not store local vector data in a contiguous
3761: array, this routine will copy the data back into the underlying
3762: vector data structure from the array obtained with `VecGetArray()`.
3764: This routine actually zeros out the `a` pointer.
3766: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3767: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3768: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3769: @*/
3770: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3771: {
3772: void *dummy;
3774: PetscFunctionBegin;
3776: PetscAssertPointer(a, 8);
3778: dummy = (void *)(*a + mstart);
3779: PetscCall(PetscFree(dummy));
3780: PetscCall(VecRestoreArrayRead(x, NULL));
3781: *a = NULL;
3782: PetscFunctionReturn(PETSC_SUCCESS);
3783: }
3785: /*@C
3786: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3787: processor's portion of the vector data. You MUST call `VecRestoreArray4dRead()`
3788: when you no longer need access to the array.
3790: Logically Collective
3792: Input Parameters:
3793: + x - the vector
3794: . m - first dimension of four dimensional array
3795: . n - second dimension of four dimensional array
3796: . p - third dimension of four dimensional array
3797: . q - fourth dimension of four dimensional array
3798: . mstart - first index you will use in first coordinate direction (often 0)
3799: . nstart - first index in the second coordinate direction (often 0)
3800: . pstart - first index in the third coordinate direction (often 0)
3801: - qstart - first index in the fourth coordinate direction (often 0)
3803: Output Parameter:
3804: . a - location to put pointer to the array
3806: Level: beginner
3808: Notes:
3809: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3810: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3811: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3812: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3814: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3816: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3817: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3818: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3819: @*/
3820: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3821: {
3822: PetscInt i, N, j, k;
3823: const PetscScalar *aa;
3824: PetscScalar ***b, **c;
3826: PetscFunctionBegin;
3828: PetscAssertPointer(a, 10);
3830: PetscCall(VecGetLocalSize(x, &N));
3831: 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);
3832: PetscCall(VecGetArrayRead(x, &aa));
3834: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3835: b = (PetscScalar ***)((*a) + m);
3836: c = (PetscScalar **)(b + m * n);
3837: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3838: for (i = 0; i < m; i++)
3839: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3840: for (i = 0; i < m; i++)
3841: for (j = 0; j < n; j++)
3842: 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;
3843: *a -= mstart;
3844: PetscFunctionReturn(PETSC_SUCCESS);
3845: }
3847: /*@C
3848: VecRestoreArray4dRead - Restores a vector after `VecGetArray4d()` has been called.
3850: Logically Collective
3852: Input Parameters:
3853: + x - the vector
3854: . m - first dimension of four dimensional array
3855: . n - second dimension of the four dimensional array
3856: . p - third dimension of the four dimensional array
3857: . q - fourth dimension of the four dimensional array
3858: . mstart - first index you will use in first coordinate direction (often 0)
3859: . nstart - first index in the second coordinate direction (often 0)
3860: . pstart - first index in the third coordinate direction (often 0)
3861: . qstart - first index in the fourth coordinate direction (often 0)
3862: - a - location of pointer to array obtained from `VecGetArray4dRead()`
3864: Level: beginner
3866: Notes:
3867: For regular PETSc vectors this routine does not involve any copies. For
3868: any special vectors that do not store local vector data in a contiguous
3869: array, this routine will copy the data back into the underlying
3870: vector data structure from the array obtained with `VecGetArray()`.
3872: This routine actually zeros out the `a` pointer.
3874: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3875: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3876: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3877: @*/
3878: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3879: {
3880: void *dummy;
3882: PetscFunctionBegin;
3884: PetscAssertPointer(a, 10);
3886: dummy = (void *)(*a + mstart);
3887: PetscCall(PetscFree(dummy));
3888: PetscCall(VecRestoreArrayRead(x, NULL));
3889: *a = NULL;
3890: PetscFunctionReturn(PETSC_SUCCESS);
3891: }
3893: #if defined(PETSC_USE_DEBUG)
3895: /*@
3896: VecLockGet - Gets the current lock status of a vector
3898: Logically Collective
3900: Input Parameter:
3901: . x - the vector
3903: Output Parameter:
3904: . state - greater than zero indicates the vector is locked for read; less than zero indicates the vector is
3905: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3907: Level: advanced
3909: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3910: @*/
3911: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3912: {
3913: PetscFunctionBegin;
3915: *state = x->lock;
3916: PetscFunctionReturn(PETSC_SUCCESS);
3917: }
3919: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3920: {
3921: PetscFunctionBegin;
3923: PetscAssertPointer(file, 2);
3924: PetscAssertPointer(func, 3);
3925: PetscAssertPointer(line, 4);
3926: #if !PetscDefined(HAVE_THREADSAFETY)
3927: {
3928: const int index = x->lockstack.currentsize - 1;
3930: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupted vec lock stack, have negative index %d", index);
3931: *file = x->lockstack.file[index];
3932: *func = x->lockstack.function[index];
3933: *line = x->lockstack.line[index];
3934: }
3935: #else
3936: *file = NULL;
3937: *func = NULL;
3938: *line = 0;
3939: #endif
3940: PetscFunctionReturn(PETSC_SUCCESS);
3941: }
3943: /*@
3944: VecLockReadPush - Pushes a read-only lock on a vector to prevent it from being written to
3946: Logically Collective
3948: Input Parameter:
3949: . x - the vector
3951: Level: intermediate
3953: Notes:
3954: If this is set then calls to `VecGetArray()` or `VecSetValues()` or any other routines that change the vectors values will generate an error.
3956: The call can be nested, i.e., called multiple times on the same vector, but each `VecLockReadPush()` has to have one matching
3957: `VecLockReadPop()`, which removes the latest read-only lock.
3959: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3960: @*/
3961: PetscErrorCode VecLockReadPush(Vec x)
3962: {
3963: PetscFunctionBegin;
3965: 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");
3966: #if !PetscDefined(HAVE_THREADSAFETY)
3967: {
3968: const char *file, *func;
3969: int index, line;
3971: if ((index = petscstack.currentsize - 2) == -1) {
3972: // vec was locked "outside" of petsc, either in user-land or main. the error message will
3973: // now show this function as the culprit, but it will include the stacktrace
3974: file = "unknown user-file";
3975: func = "unknown_user_function";
3976: line = 0;
3977: } else {
3978: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected petscstack, have negative index %d", index);
3979: file = petscstack.file[index];
3980: func = petscstack.function[index];
3981: line = petscstack.line[index];
3982: }
3983: PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
3984: }
3985: #endif
3986: PetscFunctionReturn(PETSC_SUCCESS);
3987: }
3989: /*@
3990: VecLockReadPop - Pops a read-only lock from a vector
3992: Logically Collective
3994: Input Parameter:
3995: . x - the vector
3997: Level: intermediate
3999: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
4000: @*/
4001: PetscErrorCode VecLockReadPop(Vec x)
4002: {
4003: PetscFunctionBegin;
4005: PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
4006: #if !PetscDefined(HAVE_THREADSAFETY)
4007: {
4008: const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];
4010: PetscStackPop_Private(x->lockstack, previous);
4011: }
4012: #endif
4013: PetscFunctionReturn(PETSC_SUCCESS);
4014: }
4016: /*@C
4017: VecLockWriteSet - Lock or unlock a vector for exclusive read/write access
4019: Logically Collective
4021: Input Parameters:
4022: + x - the vector
4023: - flg - `PETSC_TRUE` to lock the vector for exclusive read/write access; `PETSC_FALSE` to unlock it.
4025: Level: intermediate
4027: Notes:
4028: The function is useful in split-phase computations, which usually have a begin phase and an end phase.
4029: One can call `VecLockWriteSet`(x,`PETSC_TRUE`) in the begin phase to lock a vector for exclusive
4030: access, and call `VecLockWriteSet`(x,`PETSC_FALSE`) in the end phase to unlock the vector from exclusive
4031: access. In this way, one is ensured no other operations can access the vector in between. The code may like
4033: .vb
4034: VecGetArray(x,&xdata); // begin phase
4035: VecLockWriteSet(v,PETSC_TRUE);
4037: Other operations, which can not access x anymore (they can access xdata, of course)
4039: VecRestoreArray(x,&vdata); // end phase
4040: VecLockWriteSet(v,PETSC_FALSE);
4041: .ve
4043: The call can not be nested on the same vector, in other words, one can not call `VecLockWriteSet`(x,`PETSC_TRUE`)
4044: again before calling `VecLockWriteSet`(v,`PETSC_FALSE`).
4046: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
4047: @*/
4048: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
4049: {
4050: PetscFunctionBegin;
4052: if (flg) {
4053: 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");
4054: 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");
4055: x->lock = -1;
4056: } else {
4057: 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");
4058: x->lock = 0;
4059: }
4060: PetscFunctionReturn(PETSC_SUCCESS);
4061: }
4063: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
4064: /*@
4065: VecLockPush - Pushes a read-only lock on a vector to prevent it from being written to
4067: Level: deprecated
4069: .seealso: [](ch_vectors), `Vec`, `VecLockReadPush()`
4070: @*/
4071: PetscErrorCode VecLockPush(Vec x)
4072: {
4073: PetscFunctionBegin;
4074: PetscCall(VecLockReadPush(x));
4075: PetscFunctionReturn(PETSC_SUCCESS);
4076: }
4078: // PetscClangLinter pragma disable: -fdoc-param-list-func-parameter-documentation
4079: /*@
4080: VecLockPop - Pops a read-only lock from a vector
4082: Level: deprecated
4084: .seealso: [](ch_vectors), `Vec`, `VecLockReadPop()`
4085: @*/
4086: PetscErrorCode VecLockPop(Vec x)
4087: {
4088: PetscFunctionBegin;
4089: PetscCall(VecLockReadPop(x));
4090: PetscFunctionReturn(PETSC_SUCCESS);
4091: }
4093: #endif