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/sfimpl.h"
6: #include "petscsystypes.h"
7: #include <petsc/private/vecimpl.h>
9: PetscInt VecGetSubVectorSavedStateId = -1;
11: #if PetscDefined(USE_DEBUG)
12: // this is a no-op '0' macro in optimized builds
13: PetscErrorCode VecValidValues_Internal(Vec vec, PetscInt argnum, PetscBool begin)
14: {
15: PetscFunctionBegin;
16: if (vec->petscnative || vec->ops->getarray) {
17: PetscInt n;
18: const PetscScalar *x;
19: PetscOffloadMask mask;
21: PetscCall(VecGetOffloadMask(vec, &mask));
22: if (!PetscOffloadHost(mask)) PetscFunctionReturn(PETSC_SUCCESS);
23: PetscCall(VecGetLocalSize(vec, &n));
24: PetscCall(VecGetArrayRead(vec, &x));
25: for (PetscInt i = 0; i < n; i++) {
26: if (begin) {
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 beginning of function: Parameter number %" PetscInt_FMT, i, argnum);
28: } else {
29: 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);
30: }
31: }
32: PetscCall(VecRestoreArrayRead(vec, &x));
33: }
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
36: #endif
38: /*@
39: VecMaxPointwiseDivide - Computes the maximum of the componentwise division `max = max_i abs(x[i]/y[i])`.
41: Logically Collective
43: Input Parameters:
44: + x - the numerators
45: - y - the denominators
47: Output Parameter:
48: . max - the result
50: Level: advanced
52: Notes:
53: `x` and `y` may be the same vector
55: if a particular `y[i]` is zero, it is treated as 1 in the above formula
57: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`
58: @*/
59: PetscErrorCode VecMaxPointwiseDivide(Vec x, Vec y, PetscReal *max)
60: {
61: PetscFunctionBegin;
67: PetscCheckSameTypeAndComm(x, 1, y, 2);
68: VecCheckSameSize(x, 1, y, 2);
69: VecCheckAssembled(x);
70: VecCheckAssembled(y);
71: PetscCall(VecLockReadPush(x));
72: PetscCall(VecLockReadPush(y));
73: PetscUseTypeMethod(x, maxpointwisedivide, y, max);
74: PetscCall(VecLockReadPop(x));
75: PetscCall(VecLockReadPop(y));
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: /*@
80: VecDot - Computes the vector dot product.
82: Collective
84: Input Parameters:
85: + x - first vector
86: - y - second vector
88: Output Parameter:
89: . val - the dot product
91: Performance Issues:
92: .vb
93: per-processor memory bandwidth
94: interprocessor latency
95: work load imbalance that causes certain processes to arrive much earlier than others
96: .ve
98: Level: intermediate
100: Notes for Users of Complex Numbers:
101: For complex vectors, `VecDot()` computes
102: $ val = (x,y) = y^H x,
103: where y^H denotes the conjugate transpose of y. Note that this corresponds to the usual "mathematicians" complex
104: inner product where the SECOND argument gets the complex conjugate. Since the `BLASdot()` complex conjugates the first
105: first argument we call the `BLASdot()` with the arguments reversed.
107: Use `VecTDot()` for the indefinite form
108: $ val = (x,y) = y^T x,
109: where y^T denotes the transpose of y.
111: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDotRealPart()`
112: @*/
113: PetscErrorCode VecDot(Vec x, Vec y, PetscScalar *val)
114: {
115: PetscFunctionBegin;
121: PetscCheckSameTypeAndComm(x, 1, y, 2);
122: VecCheckSameSize(x, 1, y, 2);
123: VecCheckAssembled(x);
124: VecCheckAssembled(y);
126: PetscCall(VecLockReadPush(x));
127: PetscCall(VecLockReadPush(y));
128: PetscCall(PetscLogEventBegin(VEC_Dot, x, y, 0, 0));
129: PetscUseTypeMethod(x, dot, y, val);
130: PetscCall(PetscLogEventEnd(VEC_Dot, x, y, 0, 0));
131: PetscCall(VecLockReadPop(x));
132: PetscCall(VecLockReadPop(y));
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
136: /*@
137: VecDotRealPart - Computes the real part of the vector dot product.
139: Collective
141: Input Parameters:
142: + x - first vector
143: - y - second vector
145: Output Parameter:
146: . val - the real part of the dot product;
148: Level: intermediate
150: Performance Issues:
151: .vb
152: per-processor memory bandwidth
153: interprocessor latency
154: work load imbalance that causes certain processes to arrive much earlier than others
155: .ve
157: Notes for Users of Complex Numbers:
158: See `VecDot()` for more details on the definition of the dot product for complex numbers
160: For real numbers this returns the same value as `VecDot()`
162: 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
163: the space R^{2n} (that is a vector of 2n components with the real or imaginary part of the complex numbers for components)
165: Developer Note:
166: This is not currently optimized to compute only the real part of the dot product.
168: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`, `VecDot()`, `VecDotNorm2()`
169: @*/
170: PetscErrorCode VecDotRealPart(Vec x, Vec y, PetscReal *val)
171: {
172: PetscScalar fdot;
174: PetscFunctionBegin;
175: PetscCall(VecDot(x, y, &fdot));
176: *val = PetscRealPart(fdot);
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: /*@
181: VecNorm - Computes the vector norm.
183: Collective
185: Input Parameters:
186: + x - the vector
187: - type - the type of the norm requested
189: Output Parameter:
190: . val - the norm
192: Values of NormType:
193: + `NORM_1` - sum_i |x[i]|
194: . `NORM_2` - sqrt(sum_i |x[i]|^2)
195: . `NORM_INFINITY` - max_i |x[i]|
196: - `NORM_1_AND_2` - computes efficiently both `NORM_1` and `NORM_2` and stores them each in an output array
198: Level: intermediate
200: Notes:
201: For complex numbers `NORM_1` will return the traditional 1 norm of the 2 norm of the complex numbers; that is the 1
202: norm of the absolute values of the complex entries. In PETSc 3.6 and earlier releases it returned the 1 norm of
203: the 1 norm of the complex entries (what is returned by the BLAS routine asum()). Both are valid norms but most
204: people expect the former.
206: This routine stashes the computed norm value, repeated calls before the vector entries are changed are then rapid since the
207: precomputed value is immediately available. Certain vector operations such as `VecSet()` store the norms so the value is
208: immediately available and does not need to be explicitly computed. `VecScale()` updates any stashed norm values, thus calls after `VecScale()`
209: do not need to explicitly recompute the norm.
211: Performance Issues:
212: + per-processor memory bandwidth - limits the speed of the computation of local portion of the norm
213: . interprocessor latency - limits the accumulation of the result across ranks, .i.e. MPI_Allreduce() time
214: . number of ranks - the time for the result will grow with the log base 2 of the number of ranks sharing the vector
215: - work load imbalance - the rank with the largest number of vector entries will limit the speed up
217: .seealso: [](ch_vectors), `Vec`, `NormType`, `VecDot()`, `VecTDot()`, `VecDotBegin()`, `VecDotEnd()`, `VecNormAvailable()`,
218: `VecNormBegin()`, `VecNormEnd()`, `NormType()`
219: @*/
220: PetscErrorCode VecNorm(Vec x, NormType type, PetscReal *val)
221: {
222: PetscBool flg = PETSC_TRUE;
224: PetscFunctionBegin;
227: VecCheckAssembled(x);
231: /* Cached data? */
232: PetscCall(VecNormAvailable(x, type, &flg, val));
233: if (flg) PetscFunctionReturn(PETSC_SUCCESS);
235: PetscCall(VecLockReadPush(x));
236: PetscCall(PetscLogEventBegin(VEC_Norm, x, 0, 0, 0));
237: PetscUseTypeMethod(x, norm, type, val);
238: PetscCall(PetscLogEventEnd(VEC_Norm, x, 0, 0, 0));
239: PetscCall(VecLockReadPop(x));
241: if (type != NORM_1_AND_2) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[type], *val));
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: /*@
246: VecNormAvailable - Returns the vector norm if it is already known. That is, it has been previously computed and cached in the vector
248: Not Collective
250: Input Parameters:
251: + x - the vector
252: - 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
253: `NORM_1_AND_2`, which computes both norms and stores them
254: in a two element array.
256: Output Parameters:
257: + available - `PETSC_TRUE` if the val returned is valid
258: - val - the norm
260: Level: intermediate
262: Performance Issues:
263: .vb
264: per-processor memory bandwidth
265: interprocessor latency
266: work load imbalance that causes certain processes to arrive much earlier than others
267: .ve
269: Developer Note:
270: `PETSC_HAVE_SLOW_BLAS_NORM2` will cause a C (loop unrolled) version of the norm to be used, rather
271: than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
272: (as opposed to vendor provided) because the FORTRAN BLAS `NRM2()` routine is very slow.
274: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecTDot()`, `VecNorm()`, `VecDotBegin()`, `VecDotEnd()`,
275: `VecNormBegin()`, `VecNormEnd()`
276: @*/
277: PetscErrorCode VecNormAvailable(Vec x, NormType type, PetscBool *available, PetscReal *val)
278: {
279: PetscFunctionBegin;
285: if (type == NORM_1_AND_2) {
286: *available = PETSC_FALSE;
287: } else {
288: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[type], *val, *available));
289: }
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: /*@
294: VecNormalize - Normalizes a vector by its 2-norm.
296: Collective
298: Input Parameter:
299: . x - the vector
301: Output Parameter:
302: . val - the vector norm before normalization. May be `NULL` if the value is not needed.
304: Level: intermediate
306: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `NORM_2`, `NormType`
307: @*/
308: PetscErrorCode VecNormalize(Vec x, PetscReal *val)
309: {
310: PetscReal norm;
312: PetscFunctionBegin;
315: PetscCall(VecSetErrorIfLocked(x, 1));
317: PetscCall(PetscLogEventBegin(VEC_Normalize, x, 0, 0, 0));
318: PetscCall(VecNorm(x, NORM_2, &norm));
319: if (norm == 0.0) PetscCall(PetscInfo(x, "Vector of zero norm can not be normalized; Returning only the zero norm\n"));
320: else if (PetscIsInfOrNanReal(norm)) PetscCall(PetscInfo(x, "Vector with Inf or Nan norm can not be normalized; Returning only the norm\n"));
321: else {
322: PetscScalar s = 1.0 / norm;
323: PetscCall(VecScale(x, s));
324: }
325: PetscCall(PetscLogEventEnd(VEC_Normalize, x, 0, 0, 0));
326: if (val) *val = norm;
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: /*@C
331: VecMax - Determines the vector component with maximum real part and its location.
333: Collective
335: Input Parameter:
336: . x - the vector
338: Output Parameters:
339: + p - the index of `val` (pass `NULL` if you don't want this) in the vector
340: - val - the maximum component
342: Level: intermediate
344: Notes:
345: Returns the value `PETSC_MIN_REAL` and negative `p` if the vector is of length 0.
347: Returns the smallest index with the maximum value
349: .seealso: [](ch_vectors), `Vec`, `VecNorm()`, `VecMin()`
350: @*/
351: PetscErrorCode VecMax(Vec x, PetscInt *p, PetscReal *val)
352: {
353: PetscFunctionBegin;
356: VecCheckAssembled(x);
359: PetscCall(VecLockReadPush(x));
360: PetscCall(PetscLogEventBegin(VEC_Max, x, 0, 0, 0));
361: PetscUseTypeMethod(x, max, p, val);
362: PetscCall(PetscLogEventEnd(VEC_Max, x, 0, 0, 0));
363: PetscCall(VecLockReadPop(x));
364: PetscFunctionReturn(PETSC_SUCCESS);
365: }
367: /*@C
368: VecMin - Determines the vector component with minimum real part and its location.
370: Collective
372: Input Parameter:
373: . x - the vector
375: Output Parameters:
376: + p - the index of `val` (pass `NULL` if you don't want this location) in the vector
377: - val - the minimum component
379: Level: intermediate
381: Notes:
382: Returns the value `PETSC_MAX_REAL` and negative `p` if the vector is of length 0.
384: This returns the smallest index with the minimum value
386: .seealso: [](ch_vectors), `Vec`, `VecMax()`
387: @*/
388: PetscErrorCode VecMin(Vec x, PetscInt *p, PetscReal *val)
389: {
390: PetscFunctionBegin;
393: VecCheckAssembled(x);
396: PetscCall(VecLockReadPush(x));
397: PetscCall(PetscLogEventBegin(VEC_Min, x, 0, 0, 0));
398: PetscUseTypeMethod(x, min, p, val);
399: PetscCall(PetscLogEventEnd(VEC_Min, x, 0, 0, 0));
400: PetscCall(VecLockReadPop(x));
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: /*@
405: VecTDot - Computes an indefinite vector dot product. That is, this
406: routine does NOT use the complex conjugate.
408: Collective
410: Input Parameters:
411: + x - first vector
412: - y - second vector
414: Output Parameter:
415: . val - the dot product
417: Level: intermediate
419: Notes for Users of Complex Numbers:
420: For complex vectors, `VecTDot()` computes the indefinite form
421: $ val = (x,y) = y^T x,
422: where y^T denotes the transpose of y.
424: Use `VecDot()` for the inner product
425: $ val = (x,y) = y^H x,
426: where y^H denotes the conjugate transpose of y.
428: .seealso: [](ch_vectors), `Vec`, `VecDot()`, `VecMTDot()`
429: @*/
430: PetscErrorCode VecTDot(Vec x, Vec y, PetscScalar *val)
431: {
432: PetscFunctionBegin;
438: PetscCheckSameTypeAndComm(x, 1, y, 2);
439: VecCheckSameSize(x, 1, y, 2);
440: VecCheckAssembled(x);
441: VecCheckAssembled(y);
443: PetscCall(VecLockReadPush(x));
444: PetscCall(VecLockReadPush(y));
445: PetscCall(PetscLogEventBegin(VEC_TDot, x, y, 0, 0));
446: PetscUseTypeMethod(x, tdot, y, val);
447: PetscCall(PetscLogEventEnd(VEC_TDot, x, y, 0, 0));
448: PetscCall(VecLockReadPop(x));
449: PetscCall(VecLockReadPop(y));
450: PetscFunctionReturn(PETSC_SUCCESS);
451: }
453: /*@
454: VecScale - Scales a vector.
456: Not Collective
458: Input Parameters:
459: + x - the vector
460: - alpha - the scalar
462: Level: intermediate
464: Note:
465: For a vector with n components, `VecScale()` computes x[i] = alpha * x[i], for i=1,...,n.
467: .seealso: [](ch_vectors), `Vec`, `VecSet()`
468: @*/
469: PetscErrorCode VecScale(Vec x, PetscScalar alpha)
470: {
471: PetscReal norms[4];
472: PetscBool flgs[4];
473: PetscScalar one = 1.0;
475: PetscFunctionBegin;
478: VecCheckAssembled(x);
479: PetscCall(VecSetErrorIfLocked(x, 1));
480: if (alpha == one) PetscFunctionReturn(PETSC_SUCCESS);
482: /* get current stashed norms */
483: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
485: PetscCall(PetscLogEventBegin(VEC_Scale, x, 0, 0, 0));
486: PetscUseTypeMethod(x, scale, alpha);
487: PetscCall(PetscLogEventEnd(VEC_Scale, x, 0, 0, 0));
489: PetscCall(PetscObjectStateIncrease((PetscObject)x));
490: /* put the scaled stashed norms back into the Vec */
491: for (PetscInt i = 0; i < 4; i++) {
492: PetscReal ar = PetscAbsScalar(alpha);
493: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], ar * norms[i]));
494: }
495: PetscFunctionReturn(PETSC_SUCCESS);
496: }
498: /*@
499: VecSet - Sets all components of a vector to a single scalar value.
501: Logically Collective
503: Input Parameters:
504: + x - the vector
505: - alpha - the scalar
507: Level: beginner
509: Notes:
510: For a vector of dimension n, `VecSet()` sets x[i] = alpha, for i=1,...,n,
511: so that all vector entries then equal the identical
512: scalar value, `alpha`. Use the more general routine
513: `VecSetValues()` to set different vector entries.
515: You CANNOT call this after you have called `VecSetValues()` but before you call
516: `VecAssemblyBegin()`
518: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecSetRandom()`
519: @*/
520: PetscErrorCode VecSet(Vec x, PetscScalar alpha)
521: {
522: PetscFunctionBegin;
525: VecCheckAssembled(x);
527: PetscCall(VecSetErrorIfLocked(x, 1));
529: PetscCall(PetscLogEventBegin(VEC_Set, x, 0, 0, 0));
530: PetscUseTypeMethod(x, set, alpha);
531: PetscCall(PetscLogEventEnd(VEC_Set, x, 0, 0, 0));
532: PetscCall(PetscObjectStateIncrease((PetscObject)x));
534: /* norms can be simply set (if |alpha|*N not too large) */
536: {
537: PetscReal val = PetscAbsScalar(alpha);
538: const PetscInt N = x->map->N;
540: if (N == 0) {
541: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], 0.0l));
542: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], 0.0));
543: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], 0.0));
544: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], 0.0));
545: } else if (val > PETSC_MAX_REAL / N) {
546: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
547: } else {
548: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_1], N * val));
549: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_INFINITY], val));
550: val *= PetscSqrtReal((PetscReal)N);
551: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_2], val));
552: PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[NORM_FROBENIUS], val));
553: }
554: }
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }
558: /*@
559: VecAXPY - Computes `y = alpha x + y`.
561: Logically Collective
563: Input Parameters:
564: + alpha - the scalar
565: . x - vector scale by `alpha`
566: - y - vector accumulated into
568: Output Parameter:
569: . y - output vector
571: Level: intermediate
573: Notes:
574: This routine is optimized for alpha of 0.0, otherwise it calls the BLAS routine
575: .vb
576: VecAXPY(y,alpha,x) y = alpha x + y
577: VecAYPX(y,beta,x) y = x + beta y
578: VecAXPBY(y,alpha,beta,x) y = alpha x + beta y
579: VecWAXPY(w,alpha,x,y) w = alpha x + y
580: VecAXPBYPCZ(w,alpha,beta,gamma,x,y) z = alpha x + beta y + gamma z
581: VecMAXPY(y,nv,alpha[],x[]) y = sum alpha[i] x[i] + y
582: .ve
584: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
585: @*/
586: PetscErrorCode VecAXPY(Vec y, PetscScalar alpha, Vec x)
587: {
588: PetscFunctionBegin;
593: PetscCheckSameTypeAndComm(x, 3, y, 1);
594: VecCheckSameSize(x, 3, y, 1);
595: VecCheckAssembled(x);
596: VecCheckAssembled(y);
598: if (alpha == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
599: PetscCall(VecSetErrorIfLocked(y, 1));
600: if (x == y) {
601: PetscCall(VecScale(y, alpha + 1.0));
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
604: PetscCall(VecLockReadPush(x));
605: PetscCall(PetscLogEventBegin(VEC_AXPY, x, y, 0, 0));
606: PetscUseTypeMethod(y, axpy, alpha, x);
607: PetscCall(PetscLogEventEnd(VEC_AXPY, x, y, 0, 0));
608: PetscCall(VecLockReadPop(x));
609: PetscCall(PetscObjectStateIncrease((PetscObject)y));
610: PetscFunctionReturn(PETSC_SUCCESS);
611: }
613: /*@
614: VecAYPX - Computes `y = x + beta y`.
616: Logically Collective
618: Input Parameters:
619: + beta - the scalar
620: . x - the unscaled vector
621: - y - the vector to be scaled
623: Output Parameter:
624: . y - output vector
626: Level: intermediate
628: Developer Note:
629: The implementation is optimized for `beta` of -1.0, 0.0, and 1.0
631: .seealso: [](ch_vectors), `Vec`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
632: @*/
633: PetscErrorCode VecAYPX(Vec y, PetscScalar beta, Vec x)
634: {
635: PetscFunctionBegin;
640: PetscCheckSameTypeAndComm(x, 3, y, 1);
641: VecCheckSameSize(x, 1, y, 3);
642: VecCheckAssembled(x);
643: VecCheckAssembled(y);
645: PetscCall(VecSetErrorIfLocked(y, 1));
646: if (x == y) {
647: PetscCall(VecScale(y, beta + 1.0));
648: PetscFunctionReturn(PETSC_SUCCESS);
649: }
650: PetscCall(VecLockReadPush(x));
651: if (beta == (PetscScalar)0.0) {
652: PetscCall(VecCopy(x, y));
653: } else {
654: PetscCall(PetscLogEventBegin(VEC_AYPX, x, y, 0, 0));
655: PetscUseTypeMethod(y, aypx, beta, x);
656: PetscCall(PetscLogEventEnd(VEC_AYPX, x, y, 0, 0));
657: PetscCall(PetscObjectStateIncrease((PetscObject)y));
658: }
659: PetscCall(VecLockReadPop(x));
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
663: /*@
664: VecAXPBY - Computes `y = alpha x + beta y`.
666: Logically Collective
668: Input Parameters:
669: + alpha - first scalar
670: . beta - second scalar
671: . x - the first scaled vector
672: - y - the second scaled vector
674: Output Parameter:
675: . y - output vector
677: Level: intermediate
679: Developer Note:
680: The implementation is optimized for `alpha` and/or `beta` values of 0.0 and 1.0
682: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`
683: @*/
684: PetscErrorCode VecAXPBY(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
685: {
686: PetscFunctionBegin;
691: PetscCheckSameTypeAndComm(x, 4, y, 1);
692: VecCheckSameSize(y, 1, x, 4);
693: VecCheckAssembled(x);
694: VecCheckAssembled(y);
697: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
698: if (x == y) {
699: PetscCall(VecScale(y, alpha + beta));
700: PetscFunctionReturn(PETSC_SUCCESS);
701: }
703: PetscCall(VecSetErrorIfLocked(y, 1));
704: PetscCall(VecLockReadPush(x));
705: PetscCall(PetscLogEventBegin(VEC_AXPY, y, x, 0, 0));
706: PetscUseTypeMethod(y, axpby, alpha, beta, x);
707: PetscCall(PetscLogEventEnd(VEC_AXPY, y, x, 0, 0));
708: PetscCall(PetscObjectStateIncrease((PetscObject)y));
709: PetscCall(VecLockReadPop(x));
710: PetscFunctionReturn(PETSC_SUCCESS);
711: }
713: /*@
714: VecAXPBYPCZ - Computes `z = alpha x + beta y + gamma z`
716: Logically Collective
718: Input Parameters:
719: + alpha - first scalar
720: . beta - second scalar
721: . gamma - third scalar
722: . x - first vector
723: . y - second vector
724: - z - third vector
726: Output Parameter:
727: . z - output vector
729: Level: intermediate
731: Note:
732: `x`, `y` and `z` must be different vectors
734: Developer Note:
735: The implementation is optimized for `alpha` of 1.0 and `gamma` of 1.0 or 0.0
737: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecMAXPY()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBY()`
738: @*/
739: PetscErrorCode VecAXPBYPCZ(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
740: {
741: PetscFunctionBegin;
748: PetscCheckSameTypeAndComm(x, 5, y, 6);
749: PetscCheckSameTypeAndComm(x, 5, z, 1);
750: VecCheckSameSize(x, 1, y, 5);
751: VecCheckSameSize(x, 1, z, 6);
752: PetscCheck(x != y && x != z, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
753: PetscCheck(y != z, PetscObjectComm((PetscObject)y), PETSC_ERR_ARG_IDN, "x, y, and z must be different vectors");
754: VecCheckAssembled(x);
755: VecCheckAssembled(y);
756: VecCheckAssembled(z);
760: if (alpha == (PetscScalar)0.0 && beta == (PetscScalar)0.0 && gamma == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
762: PetscCall(VecSetErrorIfLocked(z, 1));
763: PetscCall(VecLockReadPush(x));
764: PetscCall(VecLockReadPush(y));
765: PetscCall(PetscLogEventBegin(VEC_AXPBYPCZ, x, y, z, 0));
766: PetscUseTypeMethod(z, axpbypcz, alpha, beta, gamma, x, y);
767: PetscCall(PetscLogEventEnd(VEC_AXPBYPCZ, x, y, z, 0));
768: PetscCall(PetscObjectStateIncrease((PetscObject)z));
769: PetscCall(VecLockReadPop(x));
770: PetscCall(VecLockReadPop(y));
771: PetscFunctionReturn(PETSC_SUCCESS);
772: }
774: /*@
775: VecWAXPY - Computes `w = alpha x + y`.
777: Logically Collective
779: Input Parameters:
780: + alpha - the scalar
781: . x - first vector, multiplied by `alpha`
782: - y - second vector
784: Output Parameter:
785: . w - the result
787: Level: intermediate
789: Note:
790: `w` cannot be either `x` or `y`, but `x` and `y` can be the same
792: Developer Note:
793: The implementation is optimized for alpha of -1.0, 0.0, and 1.0
795: .seealso: [](ch_vectors), `Vec`, `VecAXPY()`, `VecAYPX()`, `VecAXPBY()`, `VecMAXPY()`, `VecAXPBYPCZ()`
796: @*/
797: PetscErrorCode VecWAXPY(Vec w, PetscScalar alpha, Vec x, Vec y)
798: {
799: PetscFunctionBegin;
806: PetscCheckSameTypeAndComm(x, 3, y, 4);
807: PetscCheckSameTypeAndComm(y, 4, w, 1);
808: VecCheckSameSize(x, 3, y, 4);
809: VecCheckSameSize(x, 3, w, 1);
810: PetscCheck(w != y, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector y, suggest VecAXPY()");
811: PetscCheck(w != x, PETSC_COMM_SELF, PETSC_ERR_SUP, "Result vector w cannot be same as input vector x, suggest VecAYPX()");
812: VecCheckAssembled(x);
813: VecCheckAssembled(y);
815: PetscCall(VecSetErrorIfLocked(w, 1));
817: PetscCall(VecLockReadPush(x));
818: PetscCall(VecLockReadPush(y));
819: if (alpha == (PetscScalar)0.0) {
820: PetscCall(VecCopy(y, w));
821: } else {
822: PetscCall(PetscLogEventBegin(VEC_WAXPY, x, y, w, 0));
823: PetscUseTypeMethod(w, waxpy, alpha, x, y);
824: PetscCall(PetscLogEventEnd(VEC_WAXPY, x, y, w, 0));
825: PetscCall(PetscObjectStateIncrease((PetscObject)w));
826: }
827: PetscCall(VecLockReadPop(x));
828: PetscCall(VecLockReadPop(y));
829: PetscFunctionReturn(PETSC_SUCCESS);
830: }
832: /*@C
833: VecSetValues - Inserts or adds values into certain locations of a vector.
835: Not Collective
837: Input Parameters:
838: + x - vector to insert in
839: . ni - number of elements to add
840: . ix - indices where to add
841: . y - array of values
842: - iora - either `INSERT_VALUES` to replace the current values or `ADD_VALUES` to add values to any existing entries
844: Level: beginner
846: Notes:
847: .vb
848: `VecSetValues()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
849: .ve
851: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
852: options cannot be mixed without intervening calls to the assembly
853: routines.
855: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
856: MUST be called after all calls to `VecSetValues()` have been completed.
858: VecSetValues() uses 0-based indices in Fortran as well as in C.
860: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
861: negative indices may be passed in ix. These rows are
862: simply ignored. This allows easily inserting element load matrices
863: with homogeneous Dirchlet boundary conditions that you don't want represented
864: in the vector.
866: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesLocal()`,
867: `VecSetValue()`, `VecSetValuesBlocked()`, `InsertMode`, `INSERT_VALUES`, `ADD_VALUES`, `VecGetValues()`
868: @*/
869: PetscErrorCode VecSetValues(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
870: {
871: PetscFunctionBeginHot;
873: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
878: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
879: PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
880: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
881: PetscCall(PetscObjectStateIncrease((PetscObject)x));
882: PetscFunctionReturn(PETSC_SUCCESS);
883: }
885: /*@C
886: VecGetValues - Gets values from certain locations of a vector. Currently
887: can only get values on the same processor on which they are owned
889: Not Collective
891: Input Parameters:
892: + x - vector to get values from
893: . ni - number of elements to get
894: - ix - indices where to get them from (in global 1d numbering)
896: Output Parameter:
897: . y - array of values
899: Level: beginner
901: Notes:
902: The user provides the allocated array y; it is NOT allocated in this routine
904: `VecGetValues()` gets y[i] = x[ix[i]], for i=0,...,ni-1.
906: `VecAssemblyBegin()` and `VecAssemblyEnd()` MUST be called before calling this if `VecSetValues()` or related routine has been called
908: VecGetValues() uses 0-based indices in Fortran as well as in C.
910: If you call `VecSetOption`(x, `VEC_IGNORE_NEGATIVE_INDICES`,`PETSC_TRUE`),
911: negative indices may be passed in ix. These rows are
912: simply ignored.
914: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`
915: @*/
916: PetscErrorCode VecGetValues(Vec x, PetscInt ni, const PetscInt ix[], PetscScalar y[])
917: {
918: PetscFunctionBegin;
920: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
924: VecCheckAssembled(x);
925: PetscUseTypeMethod(x, getvalues, ni, ix, y);
926: PetscFunctionReturn(PETSC_SUCCESS);
927: }
929: /*@C
930: VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.
932: Not Collective
934: Input Parameters:
935: + x - vector to insert in
936: . ni - number of blocks to add
937: . ix - indices where to add in block count, rather than element count
938: . y - array of values
939: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES`, adds values to any existing entries
941: Level: intermediate
943: Notes:
944: `VecSetValuesBlocked()` sets x[bs*ix[i]+j] = y[bs*i+j],
945: for j=0,...,bs-1, for i=0,...,ni-1. where bs was set with VecSetBlockSize().
947: Calls to `VecSetValuesBlocked()` with the `INSERT_VALUES` and `ADD_VALUES`
948: options cannot be mixed without intervening calls to the assembly
949: routines.
951: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
952: MUST be called after all calls to `VecSetValuesBlocked()` have been completed.
954: `VecSetValuesBlocked()` uses 0-based indices in Fortran as well as in C.
956: Negative indices may be passed in ix, these rows are
957: simply ignored. This allows easily inserting element load matrices
958: with homogeneous Dirchlet boundary conditions that you don't want represented
959: in the vector.
961: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValuesBlockedLocal()`,
962: `VecSetValues()`
963: @*/
964: PetscErrorCode VecSetValuesBlocked(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
965: {
966: PetscFunctionBeginHot;
968: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
973: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
974: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
975: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
976: PetscCall(PetscObjectStateIncrease((PetscObject)x));
977: PetscFunctionReturn(PETSC_SUCCESS);
978: }
980: /*@C
981: VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
982: using a local ordering of the nodes.
984: Not Collective
986: Input Parameters:
987: + x - vector to insert in
988: . ni - number of elements to add
989: . ix - indices where to add
990: . y - array of values
991: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
993: Level: intermediate
995: Notes:
996: `VecSetValuesLocal()` sets x[ix[i]] = y[i], for i=0,...,ni-1.
998: Calls to `VecSetValues()` with the `INSERT_VALUES` and `ADD_VALUES`
999: options cannot be mixed without intervening calls to the assembly
1000: routines.
1002: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1003: MUST be called after all calls to `VecSetValuesLocal()` have been completed.
1005: `VecSetValuesLocal()` uses 0-based indices in Fortran as well as in C.
1007: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetLocalToGlobalMapping()`,
1008: `VecSetValuesBlockedLocal()`
1009: @*/
1010: PetscErrorCode VecSetValuesLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1011: {
1012: PetscInt lixp[128], *lix = lixp;
1014: PetscFunctionBeginHot;
1016: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1021: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1022: if (!x->ops->setvalueslocal) {
1023: if (x->map->mapping) {
1024: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1025: PetscCall(ISLocalToGlobalMappingApply(x->map->mapping, ni, (PetscInt *)ix, lix));
1026: PetscUseTypeMethod(x, setvalues, ni, lix, y, iora);
1027: if (ni > 128) PetscCall(PetscFree(lix));
1028: } else PetscUseTypeMethod(x, setvalues, ni, ix, y, iora);
1029: } else PetscUseTypeMethod(x, setvalueslocal, ni, ix, y, iora);
1030: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1031: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1032: PetscFunctionReturn(PETSC_SUCCESS);
1033: }
1035: /*@
1036: VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
1037: using a local ordering of the nodes.
1039: Not Collective
1041: Input Parameters:
1042: + x - vector to insert in
1043: . ni - number of blocks to add
1044: . ix - indices where to add in block count, not element count
1045: . y - array of values
1046: - iora - either `INSERT_VALUES` replaces existing entries with new values, `ADD_VALUES` adds values to any existing entries
1048: Level: intermediate
1050: Notes:
1051: `VecSetValuesBlockedLocal()` sets x[bs*ix[i]+j] = y[bs*i+j],
1052: for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with `VecSetBlockSize()`.
1054: Calls to `VecSetValuesBlockedLocal()` with the `INSERT_VALUES` and `ADD_VALUES`
1055: options cannot be mixed without intervening calls to the assembly
1056: routines.
1058: These values may be cached, so `VecAssemblyBegin()` and `VecAssemblyEnd()`
1059: MUST be called after all calls to `VecSetValuesBlockedLocal()` have been completed.
1061: `VecSetValuesBlockedLocal()` uses 0-based indices in Fortran as well as in C.
1063: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesBlocked()`,
1064: `VecSetLocalToGlobalMapping()`
1065: @*/
1066: PetscErrorCode VecSetValuesBlockedLocal(Vec x, PetscInt ni, const PetscInt ix[], const PetscScalar y[], InsertMode iora)
1067: {
1068: PetscInt lixp[128], *lix = lixp;
1070: PetscFunctionBeginHot;
1072: if (!ni) PetscFunctionReturn(PETSC_SUCCESS);
1076: PetscCall(PetscLogEventBegin(VEC_SetValues, x, 0, 0, 0));
1077: if (x->map->mapping) {
1078: if (ni > 128) PetscCall(PetscMalloc1(ni, &lix));
1079: PetscCall(ISLocalToGlobalMappingApplyBlock(x->map->mapping, ni, (PetscInt *)ix, lix));
1080: PetscUseTypeMethod(x, setvaluesblocked, ni, lix, y, iora);
1081: if (ni > 128) PetscCall(PetscFree(lix));
1082: } else {
1083: PetscUseTypeMethod(x, setvaluesblocked, ni, ix, y, iora);
1084: }
1085: PetscCall(PetscLogEventEnd(VEC_SetValues, x, 0, 0, 0));
1086: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1087: PetscFunctionReturn(PETSC_SUCCESS);
1088: }
1090: static PetscErrorCode VecMXDot_Private(Vec x, PetscInt nv, const Vec y[], PetscScalar result[], PetscErrorCode (*mxdot)(Vec, PetscInt, const Vec[], PetscScalar[]), PetscLogEvent event)
1091: {
1092: PetscFunctionBegin;
1095: VecCheckAssembled(x);
1097: if (!nv) PetscFunctionReturn(PETSC_SUCCESS);
1099: for (PetscInt i = 0; i < nv; ++i) {
1102: PetscCheckSameTypeAndComm(x, 1, y[i], 3);
1103: VecCheckSameSize(x, 1, y[i], 3);
1104: VecCheckAssembled(y[i]);
1105: PetscCall(VecLockReadPush(y[i]));
1106: }
1110: PetscCall(VecLockReadPush(x));
1111: PetscCall(PetscLogEventBegin(event, x, *y, 0, 0));
1112: PetscCall((*mxdot)(x, nv, y, result));
1113: PetscCall(PetscLogEventEnd(event, x, *y, 0, 0));
1114: PetscCall(VecLockReadPop(x));
1115: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(y[i]));
1116: PetscFunctionReturn(PETSC_SUCCESS);
1117: }
1119: /*@
1120: VecMTDot - Computes indefinite vector multiple dot products.
1121: That is, it does NOT use the complex conjugate.
1123: Collective
1125: Input Parameters:
1126: + x - one vector
1127: . nv - number of vectors
1128: - y - array of vectors. Note that vectors are pointers
1130: Output Parameter:
1131: . val - array of the dot products
1133: Level: intermediate
1135: Notes for Users of Complex Numbers:
1136: For complex vectors, `VecMTDot()` computes the indefinite form
1137: $ val = (x,y) = y^T x,
1138: where y^T denotes the transpose of y.
1140: Use `VecMDot()` for the inner product
1141: $ val = (x,y) = y^H x,
1142: where y^H denotes the conjugate transpose of y.
1144: .seealso: [](ch_vectors), `Vec`, `VecMDot()`, `VecTDot()`
1145: @*/
1146: PetscErrorCode VecMTDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1147: {
1148: PetscFunctionBegin;
1150: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mtdot, VEC_MTDot));
1151: PetscFunctionReturn(PETSC_SUCCESS);
1152: }
1154: /*@
1155: VecMDot - Computes multiple vector dot products.
1157: Collective
1159: Input Parameters:
1160: + x - one vector
1161: . nv - number of vectors
1162: - y - array of vectors.
1164: Output Parameter:
1165: . val - array of the dot products (does not allocate the array)
1167: Level: intermediate
1169: Notes for Users of Complex Numbers:
1170: For complex vectors, `VecMDot()` computes
1171: $ val = (x,y) = y^H x,
1172: where y^H denotes the conjugate transpose of y.
1174: Use `VecMTDot()` for the indefinite form
1175: $ val = (x,y) = y^T x,
1176: where y^T denotes the transpose of y.
1178: .seealso: [](ch_vectors), `Vec`, `VecMTDot()`, `VecDot()`
1179: @*/
1180: PetscErrorCode VecMDot(Vec x, PetscInt nv, const Vec y[], PetscScalar val[])
1181: {
1182: PetscFunctionBegin;
1184: PetscCall(VecMXDot_Private(x, nv, y, val, x->ops->mdot, VEC_MDot));
1185: PetscFunctionReturn(PETSC_SUCCESS);
1186: }
1188: /*@
1189: VecMAXPY - Computes `y = y + sum alpha[i] x[i]`
1191: Logically Collective
1193: Input Parameters:
1194: + nv - number of scalars and x-vectors
1195: . alpha - array of scalars
1196: . y - one vector
1197: - x - array of vectors
1199: Level: intermediate
1201: Note:
1202: `y` cannot be any of the `x` vectors
1204: .seealso: [](ch_vectors), `Vec`, `VecAYPX()`, `VecWAXPY()`, `VecAXPY()`, `VecAXPBYPCZ()`, `VecAXPBY()`
1205: @*/
1206: PetscErrorCode VecMAXPY(Vec y, PetscInt nv, const PetscScalar alpha[], Vec x[])
1207: {
1208: PetscFunctionBegin;
1210: VecCheckAssembled(y);
1212: PetscCall(VecSetErrorIfLocked(y, 1));
1213: PetscCheck(nv >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vectors (given %" PetscInt_FMT ") cannot be negative", nv);
1214: if (nv) {
1215: PetscInt zeros = 0;
1219: for (PetscInt i = 0; i < nv; ++i) {
1223: PetscCheckSameTypeAndComm(y, 1, x[i], 4);
1224: VecCheckSameSize(y, 1, x[i], 4);
1225: PetscCheck(y != x[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Array of vectors 'x' cannot contain y, found x[%" PetscInt_FMT "] == y", i);
1226: VecCheckAssembled(x[i]);
1227: PetscCall(VecLockReadPush(x[i]));
1228: zeros += alpha[i] == (PetscScalar)0.0;
1229: }
1231: if (zeros < nv) {
1232: PetscCall(PetscLogEventBegin(VEC_MAXPY, y, *x, 0, 0));
1233: PetscUseTypeMethod(y, maxpy, nv, alpha, x);
1234: PetscCall(PetscLogEventEnd(VEC_MAXPY, y, *x, 0, 0));
1235: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1236: }
1238: for (PetscInt i = 0; i < nv; ++i) PetscCall(VecLockReadPop(x[i]));
1239: }
1240: PetscFunctionReturn(PETSC_SUCCESS);
1241: }
1243: /*@
1244: VecConcatenate - Creates a new vector that is a vertical concatenation of all the given array of vectors
1245: in the order they appear in the array. The concatenated vector resides on the same
1246: communicator and is the same type as the source vectors.
1248: Collective
1250: Input Parameters:
1251: + nx - number of vectors to be concatenated
1252: - X - array containing the vectors to be concatenated in the order of concatenation
1254: Output Parameters:
1255: + Y - concatenated vector
1256: - x_is - array of index sets corresponding to the concatenated components of `Y` (pass `NULL` if not needed)
1258: Level: advanced
1260: Notes:
1261: Concatenation is similar to the functionality of a `VECNEST` object; they both represent combination of
1262: different vector spaces. However, concatenated vectors do not store any information about their
1263: sub-vectors and own their own data. Consequently, this function provides index sets to enable the
1264: manipulation of data in the concatenated vector that corresponds to the original components at creation.
1266: This is a useful tool for outer loop algorithms, particularly constrained optimizers, where the solver
1267: has to operate on combined vector spaces and cannot utilize `VECNEST` objects due to incompatibility with
1268: bound projections.
1270: .seealso: [](ch_vectors), `Vec`, `VECNEST`, `VECSCATTER`, `VecScatterCreate()`
1271: @*/
1272: PetscErrorCode VecConcatenate(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
1273: {
1274: MPI_Comm comm;
1275: VecType vec_type;
1276: Vec Ytmp, Xtmp;
1277: IS *is_tmp;
1278: PetscInt i, shift = 0, Xnl, Xng, Xbegin;
1280: PetscFunctionBegin;
1286: if ((*X)->ops->concatenate) {
1287: /* use the dedicated concatenation function if available */
1288: PetscCall((*(*X)->ops->concatenate)(nx, X, Y, x_is));
1289: } else {
1290: /* loop over vectors and start creating IS */
1291: comm = PetscObjectComm((PetscObject)(*X));
1292: PetscCall(VecGetType(*X, &vec_type));
1293: PetscCall(PetscMalloc1(nx, &is_tmp));
1294: for (i = 0; i < nx; i++) {
1295: PetscCall(VecGetSize(X[i], &Xng));
1296: PetscCall(VecGetLocalSize(X[i], &Xnl));
1297: PetscCall(VecGetOwnershipRange(X[i], &Xbegin, NULL));
1298: PetscCall(ISCreateStride(comm, Xnl, shift + Xbegin, 1, &is_tmp[i]));
1299: shift += Xng;
1300: }
1301: /* create the concatenated vector */
1302: PetscCall(VecCreate(comm, &Ytmp));
1303: PetscCall(VecSetType(Ytmp, vec_type));
1304: PetscCall(VecSetSizes(Ytmp, PETSC_DECIDE, shift));
1305: PetscCall(VecSetUp(Ytmp));
1306: /* copy data from X array to Y and return */
1307: for (i = 0; i < nx; i++) {
1308: PetscCall(VecGetSubVector(Ytmp, is_tmp[i], &Xtmp));
1309: PetscCall(VecCopy(X[i], Xtmp));
1310: PetscCall(VecRestoreSubVector(Ytmp, is_tmp[i], &Xtmp));
1311: }
1312: *Y = Ytmp;
1313: if (x_is) {
1314: *x_is = is_tmp;
1315: } else {
1316: for (i = 0; i < nx; i++) PetscCall(ISDestroy(&is_tmp[i]));
1317: PetscCall(PetscFree(is_tmp));
1318: }
1319: }
1320: PetscFunctionReturn(PETSC_SUCCESS);
1321: }
1323: /* A helper function for VecGetSubVector to check if we can implement it with no-copy (i.e. the subvector shares
1324: memory with the original vector), and the block size of the subvector.
1326: Input Parameters:
1327: + X - the original vector
1328: - is - the index set of the subvector
1330: Output Parameters:
1331: + contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
1332: . start - start of contiguous block, as an offset from the start of the ownership range of the original vector
1333: - blocksize - the block size of the subvector
1335: */
1336: PetscErrorCode VecGetSubVectorContiguityAndBS_Private(Vec X, IS is, PetscBool *contig, PetscInt *start, PetscInt *blocksize)
1337: {
1338: PetscInt gstart, gend, lstart;
1339: PetscBool red[2] = {PETSC_TRUE /*contiguous*/, PETSC_TRUE /*validVBS*/};
1340: PetscInt n, N, ibs, vbs, bs = -1;
1342: PetscFunctionBegin;
1343: PetscCall(ISGetLocalSize(is, &n));
1344: PetscCall(ISGetSize(is, &N));
1345: PetscCall(ISGetBlockSize(is, &ibs));
1346: PetscCall(VecGetBlockSize(X, &vbs));
1347: PetscCall(VecGetOwnershipRange(X, &gstart, &gend));
1348: PetscCall(ISContiguousLocal(is, gstart, gend, &lstart, &red[0]));
1349: /* block size is given by IS if ibs > 1; otherwise, check the vector */
1350: if (ibs > 1) {
1351: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1352: bs = ibs;
1353: } else {
1354: if (n % vbs || vbs == 1) red[1] = PETSC_FALSE; /* this process invalidate the collectiveness of block size */
1355: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, red, 2, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is)));
1356: if (red[0] && red[1]) bs = vbs; /* all processes have a valid block size and the access will be contiguous */
1357: }
1359: *contig = red[0];
1360: *start = lstart;
1361: *blocksize = bs;
1362: PetscFunctionReturn(PETSC_SUCCESS);
1363: }
1365: /* A helper function for VecGetSubVector, to be used when we have to build a standalone subvector through VecScatter
1367: Input Parameters:
1368: + X - the original vector
1369: . is - the index set of the subvector
1370: - bs - the block size of the subvector, gotten from VecGetSubVectorContiguityAndBS_Private()
1372: Output Parameter:
1373: . Z - the subvector, which will compose the VecScatter context on output
1374: */
1375: PetscErrorCode VecGetSubVectorThroughVecScatter_Private(Vec X, IS is, PetscInt bs, Vec *Z)
1376: {
1377: PetscInt n, N;
1378: VecScatter vscat;
1379: Vec Y;
1381: PetscFunctionBegin;
1382: PetscCall(ISGetLocalSize(is, &n));
1383: PetscCall(ISGetSize(is, &N));
1384: PetscCall(VecCreate(PetscObjectComm((PetscObject)is), &Y));
1385: PetscCall(VecSetSizes(Y, n, N));
1386: PetscCall(VecSetBlockSize(Y, bs));
1387: PetscCall(VecSetType(Y, ((PetscObject)X)->type_name));
1388: PetscCall(VecScatterCreate(X, is, Y, NULL, &vscat));
1389: PetscCall(VecScatterBegin(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1390: PetscCall(VecScatterEnd(vscat, X, Y, INSERT_VALUES, SCATTER_FORWARD));
1391: PetscCall(PetscObjectCompose((PetscObject)Y, "VecGetSubVector_Scatter", (PetscObject)vscat));
1392: PetscCall(VecScatterDestroy(&vscat));
1393: *Z = Y;
1394: PetscFunctionReturn(PETSC_SUCCESS);
1395: }
1397: /*@
1398: VecGetSubVector - Gets a vector representing part of another vector
1400: Collective
1402: Input Parameters:
1403: + X - vector from which to extract a subvector
1404: - is - index set representing portion of `X` to extract
1406: Output Parameter:
1407: . Y - subvector corresponding to `is`
1409: Level: advanced
1411: Notes:
1412: The subvector `Y` should be returned with `VecRestoreSubVector()`.
1413: `X` and must be defined on the same communicator
1415: This function may return a subvector without making a copy, therefore it is not safe to use the original vector while
1416: modifying the subvector. Other non-overlapping subvectors can still be obtained from X using this function.
1418: 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`.
1420: .seealso: [](ch_vectors), `Vec`, `IS`, `VECNEST`, `MatCreateSubMatrix()`
1421: @*/
1422: PetscErrorCode VecGetSubVector(Vec X, IS is, Vec *Y)
1423: {
1424: Vec Z;
1426: PetscFunctionBegin;
1429: PetscCheckSameComm(X, 1, is, 2);
1431: if (X->ops->getsubvector) {
1432: PetscUseTypeMethod(X, getsubvector, is, &Z);
1433: } else { /* Default implementation currently does no caching */
1434: PetscBool contig;
1435: PetscInt n, N, start, bs;
1437: PetscCall(ISGetLocalSize(is, &n));
1438: PetscCall(ISGetSize(is, &N));
1439: PetscCall(VecGetSubVectorContiguityAndBS_Private(X, is, &contig, &start, &bs));
1440: if (contig) { /* We can do a no-copy implementation */
1441: const PetscScalar *x;
1442: PetscInt state = 0;
1443: PetscBool isstd, iscuda, iship;
1445: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &isstd, VECSEQ, VECMPI, VECSTANDARD, ""));
1446: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1447: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1448: if (iscuda) {
1449: #if defined(PETSC_HAVE_CUDA)
1450: const PetscScalar *x_d;
1451: PetscMPIInt size;
1452: PetscOffloadMask flg;
1454: PetscCall(VecCUDAGetArrays_Private(X, &x, &x_d, &flg));
1455: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1456: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1457: if (x) x += start;
1458: if (x_d) x_d += start;
1459: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1460: if (size == 1) {
1461: PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1462: } else {
1463: PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1464: }
1465: Z->offloadmask = flg;
1466: #endif
1467: } else if (iship) {
1468: #if defined(PETSC_HAVE_HIP)
1469: const PetscScalar *x_d;
1470: PetscMPIInt size;
1471: PetscOffloadMask flg;
1473: PetscCall(VecHIPGetArrays_Private(X, &x, &x_d, &flg));
1474: PetscCheck(flg != PETSC_OFFLOAD_UNALLOCATED, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for PETSC_OFFLOAD_UNALLOCATED");
1475: PetscCheck(!n || x || x_d, PETSC_COMM_SELF, PETSC_ERR_SUP, "Missing vector data");
1476: if (x) x += start;
1477: if (x_d) x_d += start;
1478: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1479: if (size == 1) {
1480: PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, x, x_d, &Z));
1481: } else {
1482: PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)X), bs, n, N, x, x_d, &Z));
1483: }
1484: Z->offloadmask = flg;
1485: #endif
1486: } else if (isstd) {
1487: PetscMPIInt size;
1489: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)X), &size));
1490: PetscCall(VecGetArrayRead(X, &x));
1491: if (x) x += start;
1492: if (size == 1) {
1493: PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)X), bs, n, x, &Z));
1494: } else {
1495: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)X), bs, n, N, x, &Z));
1496: }
1497: PetscCall(VecRestoreArrayRead(X, &x));
1498: } else { /* default implementation: use place array */
1499: PetscCall(VecGetArrayRead(X, &x));
1500: PetscCall(VecCreate(PetscObjectComm((PetscObject)X), &Z));
1501: PetscCall(VecSetType(Z, ((PetscObject)X)->type_name));
1502: PetscCall(VecSetSizes(Z, n, N));
1503: PetscCall(VecSetBlockSize(Z, bs));
1504: PetscCall(VecPlaceArray(Z, x ? x + start : NULL));
1505: PetscCall(VecRestoreArrayRead(X, &x));
1506: }
1508: /* this is relevant only in debug mode */
1509: PetscCall(VecLockGet(X, &state));
1510: if (state) PetscCall(VecLockReadPush(Z));
1511: Z->ops->placearray = NULL;
1512: Z->ops->replacearray = NULL;
1513: } else { /* Have to create a scatter and do a copy */
1514: PetscCall(VecGetSubVectorThroughVecScatter_Private(X, is, bs, &Z));
1515: }
1516: }
1517: /* Record the state when the subvector was gotten so we know whether its values need to be put back */
1518: if (VecGetSubVectorSavedStateId < 0) PetscCall(PetscObjectComposedDataRegister(&VecGetSubVectorSavedStateId));
1519: PetscCall(PetscObjectComposedDataSetInt((PetscObject)Z, VecGetSubVectorSavedStateId, 1));
1520: *Y = Z;
1521: PetscFunctionReturn(PETSC_SUCCESS);
1522: }
1524: /*@
1525: VecRestoreSubVector - Restores a subvector extracted using `VecGetSubVector()`
1527: Collective
1529: Input Parameters:
1530: + X - vector from which subvector was obtained
1531: . is - index set representing the subset of `X`
1532: - Y - subvector being restored
1534: Level: advanced
1536: .seealso: [](ch_vectors), `Vec`, `IS`, `VecGetSubVector()`
1537: @*/
1538: PetscErrorCode VecRestoreSubVector(Vec X, IS is, Vec *Y)
1539: {
1540: PETSC_UNUSED PetscObjectState dummystate = 0;
1541: PetscBool unchanged;
1543: PetscFunctionBegin;
1546: PetscCheckSameComm(X, 1, is, 2);
1550: if (X->ops->restoresubvector) PetscUseTypeMethod(X, restoresubvector, is, Y);
1551: else {
1552: PetscCall(PetscObjectComposedDataGetInt((PetscObject)*Y, VecGetSubVectorSavedStateId, dummystate, unchanged));
1553: if (!unchanged) { /* If Y's state has not changed since VecGetSubVector(), we only need to destroy Y */
1554: VecScatter scatter;
1555: PetscInt state;
1557: PetscCall(VecLockGet(X, &state));
1558: PetscCheck(state == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vec X is locked for read-only or read/write access");
1560: PetscCall(PetscObjectQuery((PetscObject)*Y, "VecGetSubVector_Scatter", (PetscObject *)&scatter));
1561: if (scatter) {
1562: PetscCall(VecScatterBegin(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1563: PetscCall(VecScatterEnd(scatter, *Y, X, INSERT_VALUES, SCATTER_REVERSE));
1564: } else {
1565: PetscBool iscuda, iship;
1566: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iscuda, VECSEQCUDA, VECMPICUDA, ""));
1567: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &iship, VECSEQHIP, VECMPIHIP, ""));
1569: if (iscuda) {
1570: #if defined(PETSC_HAVE_CUDA)
1571: PetscOffloadMask ymask = (*Y)->offloadmask;
1573: /* The offloadmask of X dictates where to move memory
1574: If X GPU data is valid, then move Y data on GPU if needed
1575: Otherwise, move back to the CPU */
1576: switch (X->offloadmask) {
1577: case PETSC_OFFLOAD_BOTH:
1578: if (ymask == PETSC_OFFLOAD_CPU) {
1579: PetscCall(VecCUDAResetArray(*Y));
1580: } else if (ymask == PETSC_OFFLOAD_GPU) {
1581: X->offloadmask = PETSC_OFFLOAD_GPU;
1582: }
1583: break;
1584: case PETSC_OFFLOAD_GPU:
1585: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecCUDAResetArray(*Y));
1586: break;
1587: case PETSC_OFFLOAD_CPU:
1588: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1589: break;
1590: case PETSC_OFFLOAD_UNALLOCATED:
1591: case PETSC_OFFLOAD_KOKKOS:
1592: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1593: }
1594: #endif
1595: } else if (iship) {
1596: #if defined(PETSC_HAVE_HIP)
1597: PetscOffloadMask ymask = (*Y)->offloadmask;
1599: /* The offloadmask of X dictates where to move memory
1600: If X GPU data is valid, then move Y data on GPU if needed
1601: Otherwise, move back to the CPU */
1602: switch (X->offloadmask) {
1603: case PETSC_OFFLOAD_BOTH:
1604: if (ymask == PETSC_OFFLOAD_CPU) {
1605: PetscCall(VecHIPResetArray(*Y));
1606: } else if (ymask == PETSC_OFFLOAD_GPU) {
1607: X->offloadmask = PETSC_OFFLOAD_GPU;
1608: }
1609: break;
1610: case PETSC_OFFLOAD_GPU:
1611: if (ymask == PETSC_OFFLOAD_CPU) PetscCall(VecHIPResetArray(*Y));
1612: break;
1613: case PETSC_OFFLOAD_CPU:
1614: if (ymask == PETSC_OFFLOAD_GPU) PetscCall(VecResetArray(*Y));
1615: break;
1616: case PETSC_OFFLOAD_UNALLOCATED:
1617: case PETSC_OFFLOAD_KOKKOS:
1618: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "This should not happen");
1619: }
1620: #endif
1621: } else {
1622: /* If OpenCL vecs updated the device memory, this triggers a copy on the CPU */
1623: PetscCall(VecResetArray(*Y));
1624: }
1625: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1626: }
1627: }
1628: }
1629: PetscCall(VecDestroy(Y));
1630: PetscFunctionReturn(PETSC_SUCCESS);
1631: }
1633: /*@
1634: VecCreateLocalVector - Creates a vector object suitable for use with `VecGetLocalVector()` and friends. You must call `VecDestroy()` when the
1635: vector is no longer needed.
1637: Not Collective.
1639: Input parameter:
1640: . v - The vector for which the local vector is desired.
1642: Output parameter:
1643: . w - Upon exit this contains the local vector.
1645: Level: beginner
1647: .seealso: [](ch_vectors), `Vec`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecRestoreLocalVector()`
1648: @*/
1649: PetscErrorCode VecCreateLocalVector(Vec v, Vec *w)
1650: {
1651: PetscMPIInt size;
1653: PetscFunctionBegin;
1656: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1657: if (size == 1) PetscCall(VecDuplicate(v, w));
1658: else if (v->ops->createlocalvector) PetscUseTypeMethod(v, createlocalvector, w);
1659: else {
1660: VecType type;
1661: PetscInt n;
1663: PetscCall(VecCreate(PETSC_COMM_SELF, w));
1664: PetscCall(VecGetLocalSize(v, &n));
1665: PetscCall(VecSetSizes(*w, n, n));
1666: PetscCall(VecGetBlockSize(v, &n));
1667: PetscCall(VecSetBlockSize(*w, n));
1668: PetscCall(VecGetType(v, &type));
1669: PetscCall(VecSetType(*w, type));
1670: }
1671: PetscFunctionReturn(PETSC_SUCCESS);
1672: }
1674: /*@
1675: VecGetLocalVectorRead - Maps the local portion of a vector into a
1676: vector.
1678: Not Collective.
1680: Input parameter:
1681: . v - The vector for which the local vector is desired.
1683: Output parameter:
1684: . w - Upon exit this contains the local vector.
1686: Level: beginner
1688: Notes:
1689: You must call `VecRestoreLocalVectorRead()` when the local
1690: vector is no longer needed.
1692: This function is similar to `VecGetArrayRead()` which maps the local
1693: portion into a raw pointer. `VecGetLocalVectorRead()` is usually
1694: almost as efficient as `VecGetArrayRead()` but in certain circumstances
1695: `VecGetLocalVectorRead()` can be much more efficient than
1696: `VecGetArrayRead()`. This is because the construction of a contiguous
1697: array representing the vector data required by `VecGetArrayRead()` can
1698: be an expensive operation for certain vector types. For example, for
1699: GPU vectors `VecGetArrayRead()` requires that the data between device
1700: and host is synchronized.
1702: Unlike `VecGetLocalVector()`, this routine is not collective and
1703: preserves cached information.
1705: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1706: @*/
1707: PetscErrorCode VecGetLocalVectorRead(Vec v, Vec w)
1708: {
1709: PetscFunctionBegin;
1712: VecCheckSameLocalSize(v, 1, w, 2);
1713: if (v->ops->getlocalvectorread) {
1714: PetscUseTypeMethod(v, getlocalvectorread, w);
1715: } else {
1716: PetscScalar *a;
1718: PetscCall(VecGetArrayRead(v, (const PetscScalar **)&a));
1719: PetscCall(VecPlaceArray(w, a));
1720: }
1721: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1722: PetscCall(VecLockReadPush(v));
1723: PetscCall(VecLockReadPush(w));
1724: PetscFunctionReturn(PETSC_SUCCESS);
1725: }
1727: /*@
1728: VecRestoreLocalVectorRead - Unmaps the local portion of a vector
1729: previously mapped into a vector using `VecGetLocalVectorRead()`.
1731: Not Collective.
1733: Input parameter:
1734: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVectorRead()`.
1735: - w - The vector into which the local portion of `v` was mapped.
1737: Level: beginner
1739: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVectorRead()`, `VecGetLocalVector()`, `VecGetArrayRead()`, `VecGetArray()`
1740: @*/
1741: PetscErrorCode VecRestoreLocalVectorRead(Vec v, Vec w)
1742: {
1743: PetscFunctionBegin;
1746: if (v->ops->restorelocalvectorread) {
1747: PetscUseTypeMethod(v, restorelocalvectorread, w);
1748: } else {
1749: const PetscScalar *a;
1751: PetscCall(VecGetArrayRead(w, &a));
1752: PetscCall(VecRestoreArrayRead(v, &a));
1753: PetscCall(VecResetArray(w));
1754: }
1755: PetscCall(VecLockReadPop(v));
1756: PetscCall(VecLockReadPop(w));
1757: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1758: PetscFunctionReturn(PETSC_SUCCESS);
1759: }
1761: /*@
1762: VecGetLocalVector - Maps the local portion of a vector into a
1763: vector.
1765: Collective
1767: Input parameter:
1768: . v - The vector for which the local vector is desired.
1770: Output parameter:
1771: . w - Upon exit this contains the local vector.
1773: Level: beginner
1775: Notes:
1776: You must call `VecRestoreLocalVector()` when the local
1777: vector is no longer needed.
1779: This function is similar to `VecGetArray()` which maps the local
1780: portion into a raw pointer. `VecGetLocalVector()` is usually about as
1781: efficient as `VecGetArray()` but in certain circumstances
1782: `VecGetLocalVector()` can be much more efficient than `VecGetArray()`.
1783: This is because the construction of a contiguous array representing
1784: the vector data required by `VecGetArray()` can be an expensive
1785: operation for certain vector types. For example, for GPU vectors
1786: `VecGetArray()` requires that the data between device and host is
1787: synchronized.
1789: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecRestoreLocalVector()`, `VecGetLocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1790: @*/
1791: PetscErrorCode VecGetLocalVector(Vec v, Vec w)
1792: {
1793: PetscFunctionBegin;
1796: VecCheckSameLocalSize(v, 1, w, 2);
1797: if (v->ops->getlocalvector) {
1798: PetscUseTypeMethod(v, getlocalvector, w);
1799: } else {
1800: PetscScalar *a;
1802: PetscCall(VecGetArray(v, &a));
1803: PetscCall(VecPlaceArray(w, a));
1804: }
1805: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1806: PetscFunctionReturn(PETSC_SUCCESS);
1807: }
1809: /*@
1810: VecRestoreLocalVector - Unmaps the local portion of a vector
1811: previously mapped into a vector using `VecGetLocalVector()`.
1813: Logically Collective.
1815: Input parameter:
1816: + v - The local portion of this vector was previously mapped into `w` using `VecGetLocalVector()`.
1817: - w - The vector into which the local portion of `v` was mapped.
1819: Level: beginner
1821: .seealso: [](ch_vectors), `Vec`, `VecCreateLocalVector()`, `VecGetLocalVector()`, `VecGetLocalVectorRead()`, `VecRestoreLocalVectorRead()`, `LocalVectorRead()`, `VecGetArrayRead()`, `VecGetArray()`
1822: @*/
1823: PetscErrorCode VecRestoreLocalVector(Vec v, Vec w)
1824: {
1825: PetscFunctionBegin;
1828: if (v->ops->restorelocalvector) {
1829: PetscUseTypeMethod(v, restorelocalvector, w);
1830: } else {
1831: PetscScalar *a;
1832: PetscCall(VecGetArray(w, &a));
1833: PetscCall(VecRestoreArray(v, &a));
1834: PetscCall(VecResetArray(w));
1835: }
1836: PetscCall(PetscObjectStateIncrease((PetscObject)w));
1837: PetscCall(PetscObjectStateIncrease((PetscObject)v));
1838: PetscFunctionReturn(PETSC_SUCCESS);
1839: }
1841: /*@C
1842: VecGetArray - Returns a pointer to a contiguous array that contains this
1843: processor's portion of the vector data. For the standard PETSc
1844: vectors, `VecGetArray()` returns a pointer to the local data array and
1845: does not use any copies. If the underlying vector data is not stored
1846: in a contiguous array this routine will copy the data to a contiguous
1847: array and return a pointer to that. You MUST call `VecRestoreArray()`
1848: when you no longer need access to the array.
1850: Logically Collective
1852: Input Parameter:
1853: . x - the vector
1855: Output Parameter:
1856: . a - location to put pointer to the array
1858: Level: beginner
1860: Fortran Note:
1861: `VecGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayF90()`
1863: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
1864: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
1865: @*/
1866: PetscErrorCode VecGetArray(Vec x, PetscScalar **a)
1867: {
1868: PetscFunctionBegin;
1870: PetscCall(VecSetErrorIfLocked(x, 1));
1871: if (x->ops->getarray) { /* The if-else order matters! VECNEST, VECCUDA etc should have ops->getarray while VECCUDA etc are petscnative */
1872: PetscUseTypeMethod(x, getarray, a);
1873: } else if (x->petscnative) { /* VECSTANDARD */
1874: *a = *((PetscScalar **)x->data);
1875: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array for vector type \"%s\"", ((PetscObject)x)->type_name);
1876: PetscFunctionReturn(PETSC_SUCCESS);
1877: }
1879: /*@C
1880: VecRestoreArray - Restores a vector after `VecGetArray()` has been called and the array is no longer needed
1882: Logically Collective
1884: Input Parameters:
1885: + x - the vector
1886: - a - location of pointer to array obtained from `VecGetArray()`
1888: Level: beginner
1890: Fortran Note:
1891: `VecRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayF90()`
1893: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
1894: `VecGetArrayPair()`, `VecRestoreArrayPair()`
1895: @*/
1896: PetscErrorCode VecRestoreArray(Vec x, PetscScalar **a)
1897: {
1898: PetscFunctionBegin;
1901: if (x->ops->restorearray) {
1902: PetscUseTypeMethod(x, restorearray, a);
1903: } else PetscCheck(x->petscnative, PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot restore array for vector type \"%s\"", ((PetscObject)x)->type_name);
1904: if (a) *a = NULL;
1905: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1906: PetscFunctionReturn(PETSC_SUCCESS);
1907: }
1908: /*@C
1909: VecGetArrayRead - Get read-only pointer to contiguous array containing this processor's portion of the vector data.
1911: Not Collective
1913: Input Parameter:
1914: . x - the vector
1916: Output Parameter:
1917: . a - the array
1919: Level: beginner
1921: Notes:
1922: The array must be returned using a matching call to `VecRestoreArrayRead()`.
1924: Unlike `VecGetArray()`, preserves cached information like vector norms.
1926: Standard PETSc vectors use contiguous storage so that this routine does not perform a copy. Other vector
1927: implementations may require a copy, but such implementations should cache the contiguous representation so that
1928: only one copy is performed when this routine is called multiple times in sequence.
1930: Fortran Note:
1931: `VecGetArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayReadF90()`
1933: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
1934: @*/
1935: PetscErrorCode VecGetArrayRead(Vec x, const PetscScalar **a)
1936: {
1937: PetscFunctionBegin;
1940: if (x->ops->getarrayread) {
1941: PetscUseTypeMethod(x, getarrayread, a);
1942: } else if (x->ops->getarray) {
1943: PetscObjectState state;
1945: /* VECNEST, VECCUDA, VECKOKKOS etc */
1946: // x->ops->getarray may bump the object state, but since we know this is a read-only get
1947: // we can just undo that
1948: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
1949: PetscUseTypeMethod(x, getarray, (PetscScalar **)a);
1950: PetscCall(PetscObjectStateSet((PetscObject)x, state));
1951: } else if (x->petscnative) {
1952: /* VECSTANDARD */
1953: *a = *((PetscScalar **)x->data);
1954: } else SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot get array read for vector type \"%s\"", ((PetscObject)x)->type_name);
1955: PetscFunctionReturn(PETSC_SUCCESS);
1956: }
1958: /*@C
1959: VecRestoreArrayRead - Restore array obtained with `VecGetArrayRead()`
1961: Not Collective
1963: Input Parameters:
1964: + vec - the vector
1965: - array - the array
1967: Level: beginner
1969: Fortran Note:
1970: `VecRestoreArrayRead()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayReadF90()`
1972: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
1973: @*/
1974: PetscErrorCode VecRestoreArrayRead(Vec x, const PetscScalar **a)
1975: {
1976: PetscFunctionBegin;
1979: if (x->petscnative) { /* VECSTANDARD, VECCUDA, VECKOKKOS etc */
1980: /* nothing */
1981: } else if (x->ops->restorearrayread) { /* VECNEST */
1982: PetscUseTypeMethod(x, restorearrayread, a);
1983: } else { /* No one? */
1984: PetscObjectState state;
1986: // x->ops->restorearray may bump the object state, but since we know this is a read-restore
1987: // we can just undo that
1988: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
1989: PetscUseTypeMethod(x, restorearray, (PetscScalar **)a);
1990: PetscCall(PetscObjectStateSet((PetscObject)x, state));
1991: }
1992: if (a) *a = NULL;
1993: PetscFunctionReturn(PETSC_SUCCESS);
1994: }
1996: /*@C
1997: VecGetArrayWrite - Returns a pointer to a contiguous array that WILL contain this
1998: processor's portion of the vector data. The values in this array are NOT valid, the caller of this
1999: routine is responsible for putting values into the array; any values it does not set will be invalid
2001: Logically Collective
2003: Input Parameter:
2004: . x - the vector
2006: Output Parameter:
2007: . a - location to put pointer to the array
2009: Level: intermediate
2011: Note:
2012: The array must be returned using a matching call to `VecRestoreArrayRead()`.
2014: For vectors associated with GPUs, the host and device vectors are not synchronized before giving access. If you need correct
2015: values in the array use `VecGetArray()`
2017: Fortran Note:
2018: `VecGetArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecGetArrayWriteF90()`
2020: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteF90()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`, `VecPlaceArray()`, `VecGetArray2d()`,
2021: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArray()`, `VecRestoreArrayWrite()`
2022: @*/
2023: PetscErrorCode VecGetArrayWrite(Vec x, PetscScalar **a)
2024: {
2025: PetscFunctionBegin;
2028: PetscCall(VecSetErrorIfLocked(x, 1));
2029: if (x->ops->getarraywrite) {
2030: PetscUseTypeMethod(x, getarraywrite, a);
2031: } else {
2032: PetscCall(VecGetArray(x, a));
2033: }
2034: PetscFunctionReturn(PETSC_SUCCESS);
2035: }
2037: /*@C
2038: VecRestoreArrayWrite - Restores a vector after `VecGetArrayWrite()` has been called.
2040: Logically Collective
2042: Input Parameters:
2043: + x - the vector
2044: - a - location of pointer to array obtained from `VecGetArray()`
2046: Level: beginner
2048: Fortran Note:
2049: `VecRestoreArrayWrite()` Fortran binding is deprecated (since PETSc 3.19), use `VecRestoreArrayWriteF90()`
2051: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteF90()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`, `VecPlaceArray()`, `VecRestoreArray2d()`,
2052: `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`
2053: @*/
2054: PetscErrorCode VecRestoreArrayWrite(Vec x, PetscScalar **a)
2055: {
2056: PetscFunctionBegin;
2059: if (x->ops->restorearraywrite) {
2060: PetscUseTypeMethod(x, restorearraywrite, a);
2061: } else if (x->ops->restorearray) {
2062: PetscUseTypeMethod(x, restorearray, a);
2063: }
2064: if (a) *a = NULL;
2065: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2066: PetscFunctionReturn(PETSC_SUCCESS);
2067: }
2069: /*@C
2070: VecGetArrays - Returns a pointer to the arrays in a set of vectors
2071: that were created by a call to `VecDuplicateVecs()`.
2073: Logically Collective; No Fortran Support
2075: Input Parameters:
2076: + x - the vectors
2077: - n - the number of vectors
2079: Output Parameter:
2080: . a - location to put pointer to the array
2082: Level: intermediate
2084: Note:
2085: You MUST call `VecRestoreArrays()` when you no longer need access to the arrays.
2087: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArrays()`
2088: @*/
2089: PetscErrorCode VecGetArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2090: {
2091: PetscInt i;
2092: PetscScalar **q;
2094: PetscFunctionBegin;
2098: PetscCheck(n > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must get at least one array n = %" PetscInt_FMT, n);
2099: PetscCall(PetscMalloc1(n, &q));
2100: for (i = 0; i < n; ++i) PetscCall(VecGetArray(x[i], &q[i]));
2101: *a = q;
2102: PetscFunctionReturn(PETSC_SUCCESS);
2103: }
2105: /*@C
2106: VecRestoreArrays - Restores a group of vectors after `VecGetArrays()`
2107: has been called.
2109: Logically Collective; No Fortran Support
2111: Input Parameters:
2112: + x - the vector
2113: . n - the number of vectors
2114: - a - location of pointer to arrays obtained from `VecGetArrays()`
2116: Notes:
2117: For regular PETSc vectors this routine does not involve any copies. For
2118: any special vectors that do not store local vector data in a contiguous
2119: array, this routine will copy the data back into the underlying
2120: vector data structure from the arrays obtained with `VecGetArrays()`.
2122: Level: intermediate
2124: .seealso: [](ch_vectors), `Vec`, `VecGetArrays()`, `VecRestoreArray()`
2125: @*/
2126: PetscErrorCode VecRestoreArrays(const Vec x[], PetscInt n, PetscScalar **a[])
2127: {
2128: PetscInt i;
2129: PetscScalar **q = *a;
2131: PetscFunctionBegin;
2136: for (i = 0; i < n; ++i) PetscCall(VecRestoreArray(x[i], &q[i]));
2137: PetscCall(PetscFree(q));
2138: PetscFunctionReturn(PETSC_SUCCESS);
2139: }
2141: /*@C
2142: VecGetArrayAndMemType - Like `VecGetArray()`, but if this is a standard device vector (e.g., `VECCUDA`), the returned pointer will be a device
2143: pointer to the device memory that contains this processor's portion of the vector data. Device data is guaranteed to have the latest value.
2144: Otherwise, when this is a host vector (e.g., `VECMPI`), this routine functions the same as `VecGetArray()` and returns a host pointer.
2146: For `VECKOKKOS`, if Kokkos is configured without device (e.g., use serial or openmp), per this function, the vector works like `VECSEQ`/`VECMPI`;
2147: otherwise, it works like `VECCUDA` or `VECHIP` etc.
2149: Logically Collective; No Fortran Support
2151: Input Parameter:
2152: . x - the vector
2154: Output Parameters:
2155: + a - location to put pointer to the array
2156: - mtype - memory type of the array
2158: Level: beginner
2160: Note:
2161: Use `VecRestoreArrayAndMemType()` when the array access is no longer needed
2163: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecGetArrayReadF90()`,
2164: `VecPlaceArray()`, `VecGetArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayWrite()`, `VecRestoreArrayWrite()`
2165: @*/
2166: PetscErrorCode VecGetArrayAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2167: {
2168: PetscFunctionBegin;
2173: PetscCall(VecSetErrorIfLocked(x, 1));
2174: if (x->ops->getarrayandmemtype) {
2175: /* VECCUDA, VECKOKKOS etc */
2176: PetscUseTypeMethod(x, getarrayandmemtype, a, mtype);
2177: } else {
2178: /* VECSTANDARD, VECNEST, VECVIENNACL */
2179: PetscCall(VecGetArray(x, a));
2180: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2181: }
2182: PetscFunctionReturn(PETSC_SUCCESS);
2183: }
2185: /*@C
2186: VecRestoreArrayAndMemType - Restores a vector after `VecGetArrayAndMemType()` has been called.
2188: Logically Collective; No Fortran Support
2190: Input Parameters:
2191: + x - the vector
2192: - a - location of pointer to array obtained from `VecGetArrayAndMemType()`
2194: Level: beginner
2196: .seealso: [](ch_vectors), `Vec`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArrayRead()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecRestoreArrayReadF90()`,
2197: `VecPlaceArray()`, `VecRestoreArray2d()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2198: @*/
2199: PetscErrorCode VecRestoreArrayAndMemType(Vec x, PetscScalar **a)
2200: {
2201: PetscFunctionBegin;
2205: if (x->ops->restorearrayandmemtype) {
2206: /* VECCUDA, VECKOKKOS etc */
2207: PetscUseTypeMethod(x, restorearrayandmemtype, a);
2208: } else {
2209: /* VECNEST, VECVIENNACL */
2210: PetscCall(VecRestoreArray(x, a));
2211: } /* VECSTANDARD does nothing */
2212: if (a) *a = NULL;
2213: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2214: PetscFunctionReturn(PETSC_SUCCESS);
2215: }
2217: /*@C
2218: VecGetArrayReadAndMemType - Like `VecGetArrayRead()`, but if the input vector is a device vector, it will return a read-only device pointer.
2219: The returned pointer is guaranteed to point to up-to-date data. For host vectors, it functions as `VecGetArrayRead()`.
2221: Not Collective; No Fortran Support
2223: Input Parameter:
2224: . x - the vector
2226: Output Parameters:
2227: + a - the array
2228: - mtype - memory type of the array
2230: Level: beginner
2232: Notes:
2233: The array must be returned using a matching call to `VecRestoreArrayReadAndMemType()`.
2235: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`, `VecGetArrayAndMemType()`
2236: @*/
2237: PetscErrorCode VecGetArrayReadAndMemType(Vec x, const PetscScalar **a, PetscMemType *mtype)
2238: {
2239: PetscFunctionBegin;
2244: if (x->ops->getarrayreadandmemtype) {
2245: /* VECCUDA/VECHIP though they are also petscnative */
2246: PetscUseTypeMethod(x, getarrayreadandmemtype, a, mtype);
2247: } else if (x->ops->getarrayandmemtype) {
2248: /* VECKOKKOS */
2249: PetscObjectState state;
2251: // see VecGetArrayRead() for why
2252: PetscCall(PetscObjectStateGet((PetscObject)x, &state));
2253: PetscUseTypeMethod(x, getarrayandmemtype, (PetscScalar **)a, mtype);
2254: PetscCall(PetscObjectStateSet((PetscObject)x, state));
2255: } else {
2256: PetscCall(VecGetArrayRead(x, a));
2257: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2258: }
2259: PetscFunctionReturn(PETSC_SUCCESS);
2260: }
2262: /*@C
2263: VecRestoreArrayReadAndMemType - Restore array obtained with `VecGetArrayReadAndMemType()`
2265: Not Collective; No Fortran Support
2267: Input Parameters:
2268: + vec - the vector
2269: - array - the array
2271: Level: beginner
2273: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadAndMemType()`, `VecRestoreArrayAndMemType()`, `VecRestoreArrayWriteAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2274: @*/
2275: PetscErrorCode VecRestoreArrayReadAndMemType(Vec x, const PetscScalar **a)
2276: {
2277: PetscFunctionBegin;
2281: if (x->ops->restorearrayreadandmemtype) {
2282: /* VECCUDA/VECHIP */
2283: PetscUseTypeMethod(x, restorearrayreadandmemtype, a);
2284: } else if (!x->petscnative) {
2285: /* VECNEST */
2286: PetscCall(VecRestoreArrayRead(x, a));
2287: }
2288: if (a) *a = NULL;
2289: PetscFunctionReturn(PETSC_SUCCESS);
2290: }
2292: /*@C
2293: VecGetArrayWriteAndMemType - Like `VecGetArrayWrite()`, but if this is a device vector it will always return
2294: a device pointer to the device memory that contains this processor's portion of the vector data.
2296: Not Collective; No Fortran Support
2298: Input Parameter:
2299: . x - the vector
2301: Output Parameters:
2302: + a - the array
2303: - mtype - memory type of the array
2305: Level: beginner
2307: Note:
2308: The array must be returned using a matching call to `VecRestoreArrayWriteAndMemType()`, where it will label the device memory as most recent.
2310: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayWriteAndMemType()`, `VecGetArrayReadAndMemType()`, `VecGetArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`,
2311: @*/
2312: PetscErrorCode VecGetArrayWriteAndMemType(Vec x, PetscScalar **a, PetscMemType *mtype)
2313: {
2314: PetscFunctionBegin;
2317: PetscCall(VecSetErrorIfLocked(x, 1));
2320: if (x->ops->getarraywriteandmemtype) {
2321: /* VECCUDA, VECHIP, VECKOKKOS etc, though they are also petscnative */
2322: PetscUseTypeMethod(x, getarraywriteandmemtype, a, mtype);
2323: } else if (x->ops->getarrayandmemtype) {
2324: PetscCall(VecGetArrayAndMemType(x, a, mtype));
2325: } else {
2326: /* VECNEST, VECVIENNACL */
2327: PetscCall(VecGetArrayWrite(x, a));
2328: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2329: }
2330: PetscFunctionReturn(PETSC_SUCCESS);
2331: }
2333: /*@C
2334: VecRestoreArrayWriteAndMemType - Restore array obtained with `VecGetArrayWriteAndMemType()`
2336: Not Collective; No Fortran Support
2338: Input Parameters:
2339: + vec - the vector
2340: - array - the array
2342: Level: beginner
2344: .seealso: [](ch_vectors), `Vec`, `VecGetArrayWriteAndMemType()`, `VecRestoreArrayAndMemType()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayPair()`, `VecRestoreArrayPair()`
2345: @*/
2346: PetscErrorCode VecRestoreArrayWriteAndMemType(Vec x, PetscScalar **a)
2347: {
2348: PetscFunctionBegin;
2351: PetscCall(VecSetErrorIfLocked(x, 1));
2353: if (x->ops->restorearraywriteandmemtype) {
2354: /* VECCUDA/VECHIP */
2355: PetscMemType PETSC_UNUSED mtype; // since this function doesn't accept a memtype?
2356: PetscUseTypeMethod(x, restorearraywriteandmemtype, a, &mtype);
2357: } else if (x->ops->restorearrayandmemtype) {
2358: PetscCall(VecRestoreArrayAndMemType(x, a));
2359: } else {
2360: PetscCall(VecRestoreArray(x, a));
2361: }
2362: if (a) *a = NULL;
2363: PetscFunctionReturn(PETSC_SUCCESS);
2364: }
2366: /*@
2367: VecPlaceArray - Allows one to replace the array in a vector with an
2368: array provided by the user. This is useful to avoid copying an array
2369: into a vector.
2371: Not Collective; No Fortran Support
2373: Input Parameters:
2374: + vec - the vector
2375: - array - the array
2377: Level: developer
2379: Notes:
2380: Use `VecReplaceArray()` instead to permanently replace the array
2382: You can return to the original array with a call to `VecResetArray()`. `vec` does not take
2383: ownership of `array` in any way.
2385: The user must free `array` themselves but be careful not to
2386: do so before the vector has either been destroyed, had its original array restored with
2387: `VecResetArray()` or permanently replaced with `VecReplaceArray()`.
2389: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
2390: @*/
2391: PetscErrorCode VecPlaceArray(Vec vec, const PetscScalar array[])
2392: {
2393: PetscFunctionBegin;
2397: PetscUseTypeMethod(vec, placearray, array);
2398: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2399: PetscFunctionReturn(PETSC_SUCCESS);
2400: }
2402: /*@C
2403: VecReplaceArray - Allows one to replace the array in a vector with an
2404: array provided by the user. This is useful to avoid copying an array
2405: into a vector.
2407: Not Collective; No Fortran Support
2409: Input Parameters:
2410: + vec - the vector
2411: - array - the array
2413: Level: developer
2415: Notes:
2416: This permanently replaces the array and frees the memory associated
2417: with the old array. Use `VecPlaceArray()` to temporarily replace the array.
2419: The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
2420: freed by the user. It will be freed when the vector is destroyed.
2422: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecPlaceArray()`, `VecResetArray()`
2423: @*/
2424: PetscErrorCode VecReplaceArray(Vec vec, const PetscScalar array[])
2425: {
2426: PetscFunctionBegin;
2429: PetscUseTypeMethod(vec, replacearray, array);
2430: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
2431: PetscFunctionReturn(PETSC_SUCCESS);
2432: }
2434: /*MC
2435: VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
2436: and makes them accessible via a Fortran pointer.
2438: Synopsis:
2439: VecDuplicateVecsF90(Vec x,PetscInt n,{Vec, pointer :: y(:)},integer ierr)
2441: Collective
2443: Input Parameters:
2444: + x - a vector to mimic
2445: - n - the number of vectors to obtain
2447: Output Parameters:
2448: + y - Fortran pointer to the array of vectors
2449: - ierr - error code
2451: Example of Usage:
2452: .vb
2453: #include <petsc/finclude/petscvec.h>
2454: use petscvec
2456: Vec x
2457: Vec, pointer :: y(:)
2458: ....
2459: call VecDuplicateVecsF90(x,2,y,ierr)
2460: call VecSet(y(2),alpha,ierr)
2461: call VecSet(y(2),alpha,ierr)
2462: ....
2463: call VecDestroyVecsF90(2,y,ierr)
2464: .ve
2466: Level: beginner
2468: Note:
2469: Use `VecDestroyVecsF90()` to free the space.
2471: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecsF90()`, `VecDuplicateVecs()`
2472: M*/
2474: /*MC
2475: VecRestoreArrayF90 - Restores a vector to a usable state after a call to
2476: `VecGetArrayF90()`.
2478: Synopsis:
2479: VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2481: Logically Collective
2483: Input Parameters:
2484: + x - vector
2485: - xx_v - the Fortran pointer to the array
2487: Output Parameter:
2488: . ierr - error code
2490: Example of Usage:
2491: .vb
2492: #include <petsc/finclude/petscvec.h>
2493: use petscvec
2495: PetscScalar, pointer :: xx_v(:)
2496: ....
2497: call VecGetArrayF90(x,xx_v,ierr)
2498: xx_v(3) = a
2499: call VecRestoreArrayF90(x,xx_v,ierr)
2500: .ve
2502: Level: beginner
2504: .seealso: [](ch_vectors), `Vec`, `VecGetArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrayReadF90()`
2505: M*/
2507: /*MC
2508: VecDestroyVecsF90 - Frees a block of vectors obtained with `VecDuplicateVecsF90()`.
2510: Synopsis:
2511: VecDestroyVecsF90(PetscInt n,{Vec, pointer :: x(:)},PetscErrorCode ierr)
2513: Collective
2515: Input Parameters:
2516: + n - the number of vectors previously obtained
2517: - x - pointer to array of vector pointers
2519: Output Parameter:
2520: . ierr - error code
2522: Level: beginner
2524: .seealso: [](ch_vectors), `Vec`, `VecDestroyVecs()`, `VecDuplicateVecsF90()`
2525: M*/
2527: /*MC
2528: VecGetArrayF90 - Accesses a vector array from Fortran. For default PETSc
2529: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2530: this routine is implementation dependent. You MUST call `VecRestoreArrayF90()`
2531: when you no longer need access to the array.
2533: Synopsis:
2534: VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2536: Logically Collective
2538: Input Parameter:
2539: . x - vector
2541: Output Parameters:
2542: + xx_v - the Fortran pointer to the array
2543: - ierr - error code
2545: Example of Usage:
2546: .vb
2547: #include <petsc/finclude/petscvec.h>
2548: use petscvec
2550: PetscScalar, pointer :: xx_v(:)
2551: ....
2552: call VecGetArrayF90(x,xx_v,ierr)
2553: xx_v(3) = a
2554: call VecRestoreArrayF90(x,xx_v,ierr)
2555: .ve
2557: Level: beginner
2559: Note:
2560: If you ONLY intend to read entries from the array and not change any entries you should use `VecGetArrayReadF90()`.
2562: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayReadF90()`
2563: M*/
2565: /*MC
2566: VecGetArrayReadF90 - Accesses a read only array from Fortran. For default PETSc
2567: vectors, `VecGetArrayF90()` returns a pointer to the local data array. Otherwise,
2568: this routine is implementation dependent. You MUST call `VecRestoreArrayReadF90()`
2569: when you no longer need access to the array.
2571: Synopsis:
2572: VecGetArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2574: Logically Collective
2576: Input Parameter:
2577: . x - vector
2579: Output Parameters:
2580: + xx_v - the Fortran pointer to the array
2581: - ierr - error code
2583: Example of Usage:
2584: .vb
2585: #include <petsc/finclude/petscvec.h>
2586: use petscvec
2588: PetscScalar, pointer :: xx_v(:)
2589: ....
2590: call VecGetArrayReadF90(x,xx_v,ierr)
2591: a = xx_v(3)
2592: call VecRestoreArrayReadF90(x,xx_v,ierr)
2593: .ve
2595: Level: beginner
2597: Note:
2598: If you intend to write entries into the array you must use `VecGetArrayF90()`.
2600: .seealso: [](ch_vectors), `Vec`, `VecRestoreArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecGetArrayF90()`
2601: M*/
2603: /*MC
2604: VecRestoreArrayReadF90 - Restores a readonly vector to a usable state after a call to
2605: `VecGetArrayReadF90()`.
2607: Synopsis:
2608: VecRestoreArrayReadF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)
2610: Logically Collective
2612: Input Parameters:
2613: + x - vector
2614: - xx_v - the Fortran pointer to the array
2616: Output Parameter:
2617: . ierr - error code
2619: Example of Usage:
2620: .vb
2621: #include <petsc/finclude/petscvec.h>
2622: use petscvec
2624: PetscScalar, pointer :: xx_v(:)
2625: ....
2626: call VecGetArrayReadF90(x,xx_v,ierr)
2627: a = xx_v(3)
2628: call VecRestoreArrayReadF90(x,xx_v,ierr)
2629: .ve
2631: Level: beginner
2633: .seealso: [](ch_vectors), `Vec`, `VecGetArrayReadF90()`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecRestoreArrayRead()`, `VecRestoreArrayF90()`
2634: M*/
2636: /*@C
2637: VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
2638: processor's portion of the vector data. You MUST call `VecRestoreArray2d()`
2639: when you no longer need access to the array.
2641: Logically Collective
2643: Input Parameters:
2644: + x - the vector
2645: . m - first dimension of two dimensional array
2646: . n - second dimension of two dimensional array
2647: . mstart - first index you will use in first coordinate direction (often 0)
2648: - nstart - first index in the second coordinate direction (often 0)
2650: Output Parameter:
2651: . a - location to put pointer to the array
2653: Level: developer
2655: Notes:
2656: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2657: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2658: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2659: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2661: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2663: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2664: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2665: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2666: @*/
2667: PetscErrorCode VecGetArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2668: {
2669: PetscInt i, N;
2670: PetscScalar *aa;
2672: PetscFunctionBegin;
2676: PetscCall(VecGetLocalSize(x, &N));
2677: 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);
2678: PetscCall(VecGetArray(x, &aa));
2680: PetscCall(PetscMalloc1(m, a));
2681: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2682: *a -= mstart;
2683: PetscFunctionReturn(PETSC_SUCCESS);
2684: }
2686: /*@C
2687: VecGetArray2dWrite - Returns a pointer to a 2d contiguous array that will contain this
2688: processor's portion of the vector data. You MUST call `VecRestoreArray2dWrite()`
2689: when you no longer need access to the array.
2691: Logically Collective
2693: Input Parameters:
2694: + x - the vector
2695: . m - first dimension of two dimensional array
2696: . n - second dimension of two dimensional array
2697: . mstart - first index you will use in first coordinate direction (often 0)
2698: - nstart - first index in the second coordinate direction (often 0)
2700: Output Parameter:
2701: . a - location to put pointer to the array
2703: Level: developer
2705: Notes:
2706: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
2707: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2708: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
2709: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
2711: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2713: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2714: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2715: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2716: @*/
2717: PetscErrorCode VecGetArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2718: {
2719: PetscInt i, N;
2720: PetscScalar *aa;
2722: PetscFunctionBegin;
2726: PetscCall(VecGetLocalSize(x, &N));
2727: 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);
2728: PetscCall(VecGetArrayWrite(x, &aa));
2730: PetscCall(PetscMalloc1(m, a));
2731: for (i = 0; i < m; i++) (*a)[i] = aa + i * n - nstart;
2732: *a -= mstart;
2733: PetscFunctionReturn(PETSC_SUCCESS);
2734: }
2736: /*@C
2737: VecRestoreArray2d - Restores a vector after `VecGetArray2d()` has been called.
2739: Logically Collective
2741: Input Parameters:
2742: + x - the vector
2743: . m - first dimension of two dimensional array
2744: . n - second dimension of the two dimensional array
2745: . mstart - first index you will use in first coordinate direction (often 0)
2746: . nstart - first index in the second coordinate direction (often 0)
2747: - a - location of pointer to array obtained from `VecGetArray2d()`
2749: Level: developer
2751: Notes:
2752: For regular PETSc vectors this routine does not involve any copies. For
2753: any special vectors that do not store local vector data in a contiguous
2754: array, this routine will copy the data back into the underlying
2755: vector data structure from the array obtained with `VecGetArray()`.
2757: This routine actually zeros out the `a` pointer.
2759: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2760: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2761: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2762: @*/
2763: PetscErrorCode VecRestoreArray2d(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2764: {
2765: void *dummy;
2767: PetscFunctionBegin;
2771: dummy = (void *)(*a + mstart);
2772: PetscCall(PetscFree(dummy));
2773: PetscCall(VecRestoreArray(x, NULL));
2774: *a = NULL;
2775: PetscFunctionReturn(PETSC_SUCCESS);
2776: }
2778: /*@C
2779: VecRestoreArray2dWrite - Restores a vector after `VecGetArray2dWrite()` has been called.
2781: Logically Collective
2783: Input Parameters:
2784: + x - the vector
2785: . m - first dimension of two dimensional array
2786: . n - second dimension of the two dimensional array
2787: . mstart - first index you will use in first coordinate direction (often 0)
2788: . nstart - first index in the second coordinate direction (often 0)
2789: - a - location of pointer to array obtained from `VecGetArray2d()`
2791: Level: developer
2793: Notes:
2794: For regular PETSc vectors this routine does not involve any copies. For
2795: any special vectors that do not store local vector data in a contiguous
2796: array, this routine will copy the data back into the underlying
2797: vector data structure from the array obtained with `VecGetArray()`.
2799: This routine actually zeros out the `a` pointer.
2801: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2802: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2803: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2804: @*/
2805: PetscErrorCode VecRestoreArray2dWrite(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
2806: {
2807: void *dummy;
2809: PetscFunctionBegin;
2813: dummy = (void *)(*a + mstart);
2814: PetscCall(PetscFree(dummy));
2815: PetscCall(VecRestoreArrayWrite(x, NULL));
2816: PetscFunctionReturn(PETSC_SUCCESS);
2817: }
2819: /*@C
2820: VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
2821: processor's portion of the vector data. You MUST call `VecRestoreArray1d()`
2822: when you no longer need access to the array.
2824: Logically Collective
2826: Input Parameters:
2827: + x - the vector
2828: . m - first dimension of two dimensional array
2829: - mstart - first index you will use in first coordinate direction (often 0)
2831: Output Parameter:
2832: . a - location to put pointer to the array
2834: Level: developer
2836: Notes:
2837: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2838: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2839: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2841: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2843: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2844: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2845: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2846: @*/
2847: PetscErrorCode VecGetArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2848: {
2849: PetscInt N;
2851: PetscFunctionBegin;
2855: PetscCall(VecGetLocalSize(x, &N));
2856: 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);
2857: PetscCall(VecGetArray(x, a));
2858: *a -= mstart;
2859: PetscFunctionReturn(PETSC_SUCCESS);
2860: }
2862: /*@C
2863: VecGetArray1dWrite - Returns a pointer to a 1d contiguous array that will contain this
2864: processor's portion of the vector data. You MUST call `VecRestoreArray1dWrite()`
2865: when you no longer need access to the array.
2867: Logically Collective
2869: Input Parameters:
2870: + x - the vector
2871: . m - first dimension of two dimensional array
2872: - mstart - first index you will use in first coordinate direction (often 0)
2874: Output Parameter:
2875: . a - location to put pointer to the array
2877: Level: developer
2879: Notes:
2880: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
2881: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2882: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
2884: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
2886: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
2887: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
2888: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2889: @*/
2890: PetscErrorCode VecGetArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2891: {
2892: PetscInt N;
2894: PetscFunctionBegin;
2898: PetscCall(VecGetLocalSize(x, &N));
2899: 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);
2900: PetscCall(VecGetArrayWrite(x, a));
2901: *a -= mstart;
2902: PetscFunctionReturn(PETSC_SUCCESS);
2903: }
2905: /*@C
2906: VecRestoreArray1d - Restores a vector after `VecGetArray1d()` has been called.
2908: Logically Collective
2910: Input Parameters:
2911: + x - the vector
2912: . m - first dimension of two dimensional array
2913: . mstart - first index you will use in first coordinate direction (often 0)
2914: - a - location of pointer to array obtained from `VecGetArray1d()`
2916: Level: developer
2918: Notes:
2919: For regular PETSc vectors this routine does not involve any copies. For
2920: any special vectors that do not store local vector data in a contiguous
2921: array, this routine will copy the data back into the underlying
2922: vector data structure from the array obtained with `VecGetArray1d()`.
2924: This routine actually zeros out the `a` pointer.
2926: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2927: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2928: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2929: @*/
2930: PetscErrorCode VecRestoreArray1d(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2931: {
2932: PetscFunctionBegin;
2935: PetscCall(VecRestoreArray(x, NULL));
2936: *a = NULL;
2937: PetscFunctionReturn(PETSC_SUCCESS);
2938: }
2940: /*@C
2941: VecRestoreArray1dWrite - Restores a vector after `VecGetArray1dWrite()` has been called.
2943: Logically Collective
2945: Input Parameters:
2946: + x - the vector
2947: . m - first dimension of two dimensional array
2948: . mstart - first index you will use in first coordinate direction (often 0)
2949: - a - location of pointer to array obtained from `VecGetArray1d()`
2951: Level: developer
2953: Notes:
2954: For regular PETSc vectors this routine does not involve any copies. For
2955: any special vectors that do not store local vector data in a contiguous
2956: array, this routine will copy the data back into the underlying
2957: vector data structure from the array obtained with `VecGetArray1d()`.
2959: This routine actually zeros out the `a` pointer.
2961: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
2962: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
2963: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
2964: @*/
2965: PetscErrorCode VecRestoreArray1dWrite(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
2966: {
2967: PetscFunctionBegin;
2970: PetscCall(VecRestoreArrayWrite(x, NULL));
2971: *a = NULL;
2972: PetscFunctionReturn(PETSC_SUCCESS);
2973: }
2975: /*@C
2976: VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
2977: processor's portion of the vector data. You MUST call `VecRestoreArray3d()`
2978: when you no longer need access to the array.
2980: Logically Collective
2982: Input Parameters:
2983: + x - the vector
2984: . m - first dimension of three dimensional array
2985: . n - second dimension of three dimensional array
2986: . p - third dimension of three dimensional array
2987: . mstart - first index you will use in first coordinate direction (often 0)
2988: . nstart - first index in the second coordinate direction (often 0)
2989: - pstart - first index in the third coordinate direction (often 0)
2991: Output Parameter:
2992: . a - location to put pointer to the array
2994: Level: developer
2996: Notes:
2997: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
2998: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
2999: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3000: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3002: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3004: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3005: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3006: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3007: @*/
3008: PetscErrorCode VecGetArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3009: {
3010: PetscInt i, N, j;
3011: PetscScalar *aa, **b;
3013: PetscFunctionBegin;
3017: PetscCall(VecGetLocalSize(x, &N));
3018: 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);
3019: PetscCall(VecGetArray(x, &aa));
3021: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3022: b = (PetscScalar **)((*a) + m);
3023: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3024: for (i = 0; i < m; i++)
3025: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3026: *a -= mstart;
3027: PetscFunctionReturn(PETSC_SUCCESS);
3028: }
3030: /*@C
3031: VecGetArray3dWrite - Returns a pointer to a 3d contiguous array that will contain this
3032: processor's portion of the vector data. You MUST call `VecRestoreArray3dWrite()`
3033: when you no longer need access to the array.
3035: Logically Collective
3037: Input Parameters:
3038: + x - the vector
3039: . m - first dimension of three dimensional array
3040: . n - second dimension of three dimensional array
3041: . p - third dimension of three dimensional array
3042: . mstart - first index you will use in first coordinate direction (often 0)
3043: . nstart - first index in the second coordinate direction (often 0)
3044: - pstart - first index in the third coordinate direction (often 0)
3046: Output Parameter:
3047: . a - location to put pointer to the array
3049: Level: developer
3051: Notes:
3052: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3053: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3054: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3055: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3057: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3059: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3060: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3061: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3062: @*/
3063: PetscErrorCode VecGetArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3064: {
3065: PetscInt i, N, j;
3066: PetscScalar *aa, **b;
3068: PetscFunctionBegin;
3072: PetscCall(VecGetLocalSize(x, &N));
3073: 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);
3074: PetscCall(VecGetArrayWrite(x, &aa));
3076: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3077: b = (PetscScalar **)((*a) + m);
3078: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3079: for (i = 0; i < m; i++)
3080: for (j = 0; j < n; j++) b[i * n + j] = aa + i * n * p + j * p - pstart;
3082: *a -= mstart;
3083: PetscFunctionReturn(PETSC_SUCCESS);
3084: }
3086: /*@C
3087: VecRestoreArray3d - Restores a vector after `VecGetArray3d()` has been called.
3089: Logically Collective
3091: Input Parameters:
3092: + x - the vector
3093: . m - first dimension of three dimensional array
3094: . n - second dimension of the three dimensional array
3095: . p - third dimension of the three dimensional array
3096: . mstart - first index you will use in first coordinate direction (often 0)
3097: . nstart - first index in the second coordinate direction (often 0)
3098: . pstart - first index in the third coordinate direction (often 0)
3099: - a - location of pointer to array obtained from VecGetArray3d()
3101: Level: developer
3103: Notes:
3104: For regular PETSc vectors this routine does not involve any copies. For
3105: any special vectors that do not store local vector data in a contiguous
3106: array, this routine will copy the data back into the underlying
3107: vector data structure from the array obtained with `VecGetArray()`.
3109: This routine actually zeros out the `a` pointer.
3111: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3112: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3113: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3114: @*/
3115: PetscErrorCode VecRestoreArray3d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3116: {
3117: void *dummy;
3119: PetscFunctionBegin;
3123: dummy = (void *)(*a + mstart);
3124: PetscCall(PetscFree(dummy));
3125: PetscCall(VecRestoreArray(x, NULL));
3126: *a = NULL;
3127: PetscFunctionReturn(PETSC_SUCCESS);
3128: }
3130: /*@C
3131: VecRestoreArray3dWrite - Restores a vector after `VecGetArray3dWrite()` has been called.
3133: Logically Collective
3135: Input Parameters:
3136: + x - the vector
3137: . m - first dimension of three dimensional array
3138: . n - second dimension of the three dimensional array
3139: . p - third dimension of the three dimensional array
3140: . mstart - first index you will use in first coordinate direction (often 0)
3141: . nstart - first index in the second coordinate direction (often 0)
3142: . pstart - first index in the third coordinate direction (often 0)
3143: - a - location of pointer to array obtained from VecGetArray3d()
3145: Level: developer
3147: Notes:
3148: For regular PETSc vectors this routine does not involve any copies. For
3149: any special vectors that do not store local vector data in a contiguous
3150: array, this routine will copy the data back into the underlying
3151: vector data structure from the array obtained with `VecGetArray()`.
3153: This routine actually zeros out the `a` pointer.
3155: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3156: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3157: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3158: @*/
3159: PetscErrorCode VecRestoreArray3dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3160: {
3161: void *dummy;
3163: PetscFunctionBegin;
3167: dummy = (void *)(*a + mstart);
3168: PetscCall(PetscFree(dummy));
3169: PetscCall(VecRestoreArrayWrite(x, NULL));
3170: *a = NULL;
3171: PetscFunctionReturn(PETSC_SUCCESS);
3172: }
3174: /*@C
3175: VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
3176: processor's portion of the vector data. You MUST call `VecRestoreArray4d()`
3177: when you no longer need access to the array.
3179: Logically Collective
3181: Input Parameters:
3182: + x - the vector
3183: . m - first dimension of four dimensional array
3184: . n - second dimension of four dimensional array
3185: . p - third dimension of four dimensional array
3186: . q - fourth dimension of four dimensional array
3187: . mstart - first index you will use in first coordinate direction (often 0)
3188: . nstart - first index in the second coordinate direction (often 0)
3189: . pstart - first index in the third coordinate direction (often 0)
3190: - qstart - first index in the fourth coordinate direction (often 0)
3192: Output Parameter:
3193: . a - location to put pointer to the array
3195: Level: developer
3197: Notes:
3198: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3199: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3200: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3201: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3203: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3205: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3206: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3207: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3208: @*/
3209: PetscErrorCode VecGetArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3210: {
3211: PetscInt i, N, j, k;
3212: PetscScalar *aa, ***b, **c;
3214: PetscFunctionBegin;
3218: PetscCall(VecGetLocalSize(x, &N));
3219: 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);
3220: PetscCall(VecGetArray(x, &aa));
3222: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3223: b = (PetscScalar ***)((*a) + m);
3224: c = (PetscScalar **)(b + m * n);
3225: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3226: for (i = 0; i < m; i++)
3227: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3228: for (i = 0; i < m; i++)
3229: for (j = 0; j < n; j++)
3230: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3231: *a -= mstart;
3232: PetscFunctionReturn(PETSC_SUCCESS);
3233: }
3235: /*@C
3236: VecGetArray4dWrite - Returns a pointer to a 4d contiguous array that will contain this
3237: processor's portion of the vector data. You MUST call `VecRestoreArray4dWrite()`
3238: when you no longer need access to the array.
3240: Logically Collective
3242: Input Parameters:
3243: + x - the vector
3244: . m - first dimension of four dimensional array
3245: . n - second dimension of four dimensional array
3246: . p - third dimension of four dimensional array
3247: . q - fourth dimension of four dimensional array
3248: . mstart - first index you will use in first coordinate direction (often 0)
3249: . nstart - first index in the second coordinate direction (often 0)
3250: . pstart - first index in the third coordinate direction (often 0)
3251: - qstart - first index in the fourth coordinate direction (often 0)
3253: Output Parameter:
3254: . a - location to put pointer to the array
3256: Level: developer
3258: Notes:
3259: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3260: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3261: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3262: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3264: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3266: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3267: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3268: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3269: @*/
3270: PetscErrorCode VecGetArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3271: {
3272: PetscInt i, N, j, k;
3273: PetscScalar *aa, ***b, **c;
3275: PetscFunctionBegin;
3279: PetscCall(VecGetLocalSize(x, &N));
3280: 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);
3281: PetscCall(VecGetArrayWrite(x, &aa));
3283: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3284: b = (PetscScalar ***)((*a) + m);
3285: c = (PetscScalar **)(b + m * n);
3286: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3287: for (i = 0; i < m; i++)
3288: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3289: for (i = 0; i < m; i++)
3290: for (j = 0; j < n; j++)
3291: for (k = 0; k < p; k++) c[i * n * p + j * p + k] = aa + i * n * p * q + j * p * q + k * q - qstart;
3292: *a -= mstart;
3293: PetscFunctionReturn(PETSC_SUCCESS);
3294: }
3296: /*@C
3297: VecRestoreArray4d - Restores a vector after `VecGetArray4d()` has been called.
3299: Logically Collective
3301: Input Parameters:
3302: + x - the vector
3303: . m - first dimension of four dimensional array
3304: . n - second dimension of the four dimensional array
3305: . p - third dimension of the four dimensional array
3306: . q - fourth dimension of the four dimensional array
3307: . mstart - first index you will use in first coordinate direction (often 0)
3308: . nstart - first index in the second coordinate direction (often 0)
3309: . pstart - first index in the third coordinate direction (often 0)
3310: . qstart - first index in the fourth coordinate direction (often 0)
3311: - a - location of pointer to array obtained from VecGetArray4d()
3313: Level: developer
3315: Notes:
3316: For regular PETSc vectors this routine does not involve any copies. For
3317: any special vectors that do not store local vector data in a contiguous
3318: array, this routine will copy the data back into the underlying
3319: vector data structure from the array obtained with `VecGetArray()`.
3321: This routine actually zeros out the `a` pointer.
3323: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3324: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3325: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3326: @*/
3327: PetscErrorCode VecRestoreArray4d(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3328: {
3329: void *dummy;
3331: PetscFunctionBegin;
3335: dummy = (void *)(*a + mstart);
3336: PetscCall(PetscFree(dummy));
3337: PetscCall(VecRestoreArray(x, NULL));
3338: *a = NULL;
3339: PetscFunctionReturn(PETSC_SUCCESS);
3340: }
3342: /*@C
3343: VecRestoreArray4dWrite - Restores a vector after `VecGetArray4dWrite()` has been called.
3345: Logically Collective
3347: Input Parameters:
3348: + x - the vector
3349: . m - first dimension of four dimensional array
3350: . n - second dimension of the four dimensional array
3351: . p - third dimension of the four dimensional array
3352: . q - fourth dimension of the four dimensional array
3353: . mstart - first index you will use in first coordinate direction (often 0)
3354: . nstart - first index in the second coordinate direction (often 0)
3355: . pstart - first index in the third coordinate direction (often 0)
3356: . qstart - first index in the fourth coordinate direction (often 0)
3357: - a - location of pointer to array obtained from `VecGetArray4d()`
3359: Level: developer
3361: Notes:
3362: For regular PETSc vectors this routine does not involve any copies. For
3363: any special vectors that do not store local vector data in a contiguous
3364: array, this routine will copy the data back into the underlying
3365: vector data structure from the array obtained with `VecGetArray()`.
3367: This routine actually zeros out the `a` pointer.
3369: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3370: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3371: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3372: @*/
3373: PetscErrorCode VecRestoreArray4dWrite(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3374: {
3375: void *dummy;
3377: PetscFunctionBegin;
3381: dummy = (void *)(*a + mstart);
3382: PetscCall(PetscFree(dummy));
3383: PetscCall(VecRestoreArrayWrite(x, NULL));
3384: *a = NULL;
3385: PetscFunctionReturn(PETSC_SUCCESS);
3386: }
3388: /*@C
3389: VecGetArray2dRead - Returns a pointer to a 2d contiguous array that contains this
3390: processor's portion of the vector data. You MUST call `VecRestoreArray2dRead()`
3391: when you no longer need access to the array.
3393: Logically Collective
3395: Input Parameters:
3396: + x - the vector
3397: . m - first dimension of two dimensional array
3398: . n - second dimension of two dimensional array
3399: . mstart - first index you will use in first coordinate direction (often 0)
3400: - nstart - first index in the second coordinate direction (often 0)
3402: Output Parameter:
3403: . a - location to put pointer to the array
3405: Level: developer
3407: Notes:
3408: For a vector obtained from `DMCreateLocalVector()` `mstart` and `nstart` are likely
3409: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3410: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3411: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray2d()`.
3413: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3415: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3416: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3417: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3418: @*/
3419: PetscErrorCode VecGetArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3420: {
3421: PetscInt i, N;
3422: const PetscScalar *aa;
3424: PetscFunctionBegin;
3428: PetscCall(VecGetLocalSize(x, &N));
3429: 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);
3430: PetscCall(VecGetArrayRead(x, &aa));
3432: PetscCall(PetscMalloc1(m, a));
3433: for (i = 0; i < m; i++) (*a)[i] = (PetscScalar *)aa + i * n - nstart;
3434: *a -= mstart;
3435: PetscFunctionReturn(PETSC_SUCCESS);
3436: }
3438: /*@C
3439: VecRestoreArray2dRead - Restores a vector after `VecGetArray2dRead()` has been called.
3441: Logically Collective
3443: Input Parameters:
3444: + x - the vector
3445: . m - first dimension of two dimensional array
3446: . n - second dimension of the two dimensional array
3447: . mstart - first index you will use in first coordinate direction (often 0)
3448: . nstart - first index in the second coordinate direction (often 0)
3449: - a - location of pointer to array obtained from VecGetArray2d()
3451: Level: developer
3453: Notes:
3454: For regular PETSc vectors this routine does not involve any copies. For
3455: any special vectors that do not store local vector data in a contiguous
3456: array, this routine will copy the data back into the underlying
3457: vector data structure from the array obtained with `VecGetArray()`.
3459: This routine actually zeros out the `a` pointer.
3461: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3462: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3463: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3464: @*/
3465: PetscErrorCode VecRestoreArray2dRead(Vec x, PetscInt m, PetscInt n, PetscInt mstart, PetscInt nstart, PetscScalar **a[])
3466: {
3467: void *dummy;
3469: PetscFunctionBegin;
3473: dummy = (void *)(*a + mstart);
3474: PetscCall(PetscFree(dummy));
3475: PetscCall(VecRestoreArrayRead(x, NULL));
3476: *a = NULL;
3477: PetscFunctionReturn(PETSC_SUCCESS);
3478: }
3480: /*@C
3481: VecGetArray1dRead - Returns a pointer to a 1d contiguous array that contains this
3482: processor's portion of the vector data. You MUST call `VecRestoreArray1dRead()`
3483: when you no longer need access to the array.
3485: Logically Collective
3487: Input Parameters:
3488: + x - the vector
3489: . m - first dimension of two dimensional array
3490: - mstart - first index you will use in first coordinate direction (often 0)
3492: Output Parameter:
3493: . a - location to put pointer to the array
3495: Level: developer
3497: Notes:
3498: For a vector obtained from `DMCreateLocalVector()` `mstart` is likely
3499: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3500: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`.
3502: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3504: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3505: `VecRestoreArray2d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3506: `VecGetArray2d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3507: @*/
3508: PetscErrorCode VecGetArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3509: {
3510: PetscInt N;
3512: PetscFunctionBegin;
3516: PetscCall(VecGetLocalSize(x, &N));
3517: 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);
3518: PetscCall(VecGetArrayRead(x, (const PetscScalar **)a));
3519: *a -= mstart;
3520: PetscFunctionReturn(PETSC_SUCCESS);
3521: }
3523: /*@C
3524: VecRestoreArray1dRead - Restores a vector after `VecGetArray1dRead()` has been called.
3526: Logically Collective
3528: Input Parameters:
3529: + x - the vector
3530: . m - first dimension of two dimensional array
3531: . mstart - first index you will use in first coordinate direction (often 0)
3532: - a - location of pointer to array obtained from `VecGetArray1dRead()`
3534: Level: developer
3536: Notes:
3537: For regular PETSc vectors this routine does not involve any copies. For
3538: any special vectors that do not store local vector data in a contiguous
3539: array, this routine will copy the data back into the underlying
3540: vector data structure from the array obtained with `VecGetArray1dRead()`.
3542: This routine actually zeros out the `a` pointer.
3544: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3545: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3546: `VecGetArray1d()`, `VecRestoreArray2d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3547: @*/
3548: PetscErrorCode VecRestoreArray1dRead(Vec x, PetscInt m, PetscInt mstart, PetscScalar *a[])
3549: {
3550: PetscFunctionBegin;
3553: PetscCall(VecRestoreArrayRead(x, NULL));
3554: *a = NULL;
3555: PetscFunctionReturn(PETSC_SUCCESS);
3556: }
3558: /*@C
3559: VecGetArray3dRead - Returns a pointer to a 3d contiguous array that contains this
3560: processor's portion of the vector data. You MUST call `VecRestoreArray3dRead()`
3561: when you no longer need access to the array.
3563: Logically Collective
3565: Input Parameters:
3566: + x - the vector
3567: . m - first dimension of three dimensional array
3568: . n - second dimension of three dimensional array
3569: . p - third dimension of three dimensional array
3570: . mstart - first index you will use in first coordinate direction (often 0)
3571: . nstart - first index in the second coordinate direction (often 0)
3572: - pstart - first index in the third coordinate direction (often 0)
3574: Output Parameter:
3575: . a - location to put pointer to the array
3577: Level: developer
3579: Notes:
3580: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3581: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3582: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3583: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3dRead()`.
3585: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3587: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3588: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3589: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3590: @*/
3591: PetscErrorCode VecGetArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3592: {
3593: PetscInt i, N, j;
3594: const PetscScalar *aa;
3595: PetscScalar **b;
3597: PetscFunctionBegin;
3601: PetscCall(VecGetLocalSize(x, &N));
3602: 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);
3603: PetscCall(VecGetArrayRead(x, &aa));
3605: PetscCall(PetscMalloc(m * sizeof(PetscScalar **) + m * n * sizeof(PetscScalar *), a));
3606: b = (PetscScalar **)((*a) + m);
3607: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3608: for (i = 0; i < m; i++)
3609: for (j = 0; j < n; j++) b[i * n + j] = (PetscScalar *)aa + i * n * p + j * p - pstart;
3610: *a -= mstart;
3611: PetscFunctionReturn(PETSC_SUCCESS);
3612: }
3614: /*@C
3615: VecRestoreArray3dRead - Restores a vector after `VecGetArray3dRead()` has been called.
3617: Logically Collective
3619: Input Parameters:
3620: + x - the vector
3621: . m - first dimension of three dimensional array
3622: . n - second dimension of the three dimensional array
3623: . p - third dimension of the three dimensional array
3624: . mstart - first index you will use in first coordinate direction (often 0)
3625: . nstart - first index in the second coordinate direction (often 0)
3626: . pstart - first index in the third coordinate direction (often 0)
3627: - a - location of pointer to array obtained from `VecGetArray3dRead()`
3629: Level: developer
3631: Notes:
3632: For regular PETSc vectors this routine does not involve any copies. For
3633: any special vectors that do not store local vector data in a contiguous
3634: array, this routine will copy the data back into the underlying
3635: vector data structure from the array obtained with `VecGetArray()`.
3637: This routine actually zeros out the `a` pointer.
3639: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3640: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3641: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3642: @*/
3643: PetscErrorCode VecRestoreArray3dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscScalar ***a[])
3644: {
3645: void *dummy;
3647: PetscFunctionBegin;
3651: dummy = (void *)(*a + mstart);
3652: PetscCall(PetscFree(dummy));
3653: PetscCall(VecRestoreArrayRead(x, NULL));
3654: *a = NULL;
3655: PetscFunctionReturn(PETSC_SUCCESS);
3656: }
3658: /*@C
3659: VecGetArray4dRead - Returns a pointer to a 4d contiguous array that contains this
3660: processor's portion of the vector data. You MUST call `VecRestoreArray4dRead()`
3661: when you no longer need access to the array.
3663: Logically Collective
3665: Input Parameters:
3666: + x - the vector
3667: . m - first dimension of four dimensional array
3668: . n - second dimension of four dimensional array
3669: . p - third dimension of four dimensional array
3670: . q - fourth dimension of four dimensional array
3671: . mstart - first index you will use in first coordinate direction (often 0)
3672: . nstart - first index in the second coordinate direction (often 0)
3673: . pstart - first index in the third coordinate direction (often 0)
3674: - qstart - first index in the fourth coordinate direction (often 0)
3676: Output Parameter:
3677: . a - location to put pointer to the array
3679: Level: beginner
3681: Notes:
3682: For a vector obtained from `DMCreateLocalVector()` `mstart`, `nstart`, and `pstart` are likely
3683: obtained from the corner indices obtained from `DMDAGetGhostCorners()` while for
3684: `DMCreateGlobalVector()` they are the corner indices from `DMDAGetCorners()`. In both cases
3685: the arguments from `DMDAGet[Ghost]Corners()` are reversed in the call to `VecGetArray3d()`.
3687: For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.
3689: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecGetArrays()`, `VecGetArrayF90()`, `VecPlaceArray()`,
3690: `VecRestoreArray2d()`, `DMDAVecGetarray()`, `DMDAVecRestoreArray()`, `VecGetArray3d()`, `VecRestoreArray3d()`,
3691: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`
3692: @*/
3693: PetscErrorCode VecGetArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3694: {
3695: PetscInt i, N, j, k;
3696: const PetscScalar *aa;
3697: PetscScalar ***b, **c;
3699: PetscFunctionBegin;
3703: PetscCall(VecGetLocalSize(x, &N));
3704: 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);
3705: PetscCall(VecGetArrayRead(x, &aa));
3707: PetscCall(PetscMalloc(m * sizeof(PetscScalar ***) + m * n * sizeof(PetscScalar **) + m * n * p * sizeof(PetscScalar *), a));
3708: b = (PetscScalar ***)((*a) + m);
3709: c = (PetscScalar **)(b + m * n);
3710: for (i = 0; i < m; i++) (*a)[i] = b + i * n - nstart;
3711: for (i = 0; i < m; i++)
3712: for (j = 0; j < n; j++) b[i * n + j] = c + i * n * p + j * p - pstart;
3713: for (i = 0; i < m; i++)
3714: for (j = 0; j < n; j++)
3715: 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;
3716: *a -= mstart;
3717: PetscFunctionReturn(PETSC_SUCCESS);
3718: }
3720: /*@C
3721: VecRestoreArray4dRead - Restores a vector after `VecGetArray4d()` has been called.
3723: Logically Collective
3725: Input Parameters:
3726: + x - the vector
3727: . m - first dimension of four dimensional array
3728: . n - second dimension of the four dimensional array
3729: . p - third dimension of the four dimensional array
3730: . q - fourth dimension of the four dimensional array
3731: . mstart - first index you will use in first coordinate direction (often 0)
3732: . nstart - first index in the second coordinate direction (often 0)
3733: . pstart - first index in the third coordinate direction (often 0)
3734: . qstart - first index in the fourth coordinate direction (often 0)
3735: - a - location of pointer to array obtained from `VecGetArray4dRead()`
3737: Level: beginner
3739: Notes:
3740: For regular PETSc vectors this routine does not involve any copies. For
3741: any special vectors that do not store local vector data in a contiguous
3742: array, this routine will copy the data back into the underlying
3743: vector data structure from the array obtained with `VecGetArray()`.
3745: This routine actually zeros out the `a` pointer.
3747: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecRestoreArrays()`, `VecRestoreArrayF90()`, `VecPlaceArray()`,
3748: `VecGetArray2d()`, `VecGetArray3d()`, `VecRestoreArray3d()`, `DMDAVecGetArray()`, `DMDAVecRestoreArray()`
3749: `VecGetArray1d()`, `VecRestoreArray1d()`, `VecGetArray4d()`, `VecRestoreArray4d()`, `VecGet`
3750: @*/
3751: PetscErrorCode VecRestoreArray4dRead(Vec x, PetscInt m, PetscInt n, PetscInt p, PetscInt q, PetscInt mstart, PetscInt nstart, PetscInt pstart, PetscInt qstart, PetscScalar ****a[])
3752: {
3753: void *dummy;
3755: PetscFunctionBegin;
3759: dummy = (void *)(*a + mstart);
3760: PetscCall(PetscFree(dummy));
3761: PetscCall(VecRestoreArrayRead(x, NULL));
3762: *a = NULL;
3763: PetscFunctionReturn(PETSC_SUCCESS);
3764: }
3766: #if defined(PETSC_USE_DEBUG)
3768: /*@
3769: VecLockGet - Gets the current lock status of a vector
3771: Logically Collective
3773: Input Parameter:
3774: . x - the vector
3776: Output Parameter:
3777: . state - greater than zero indicates the vector is locked for read; less then zero indicates the vector is
3778: locked for write; equal to zero means the vector is unlocked, that is, it is free to read or write.
3780: Level: advanced
3782: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`
3783: @*/
3784: PetscErrorCode VecLockGet(Vec x, PetscInt *state)
3785: {
3786: PetscFunctionBegin;
3788: *state = x->lock;
3789: PetscFunctionReturn(PETSC_SUCCESS);
3790: }
3792: PetscErrorCode VecLockGetLocation(Vec x, const char *file[], const char *func[], int *line)
3793: {
3794: PetscFunctionBegin;
3799: #if !PetscDefined(HAVE_THREADSAFETY)
3800: {
3801: const int index = x->lockstack.currentsize - 1;
3803: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupted vec lock stack, have negative index %d", index);
3804: *file = x->lockstack.file[index];
3805: *func = x->lockstack.function[index];
3806: *line = x->lockstack.line[index];
3807: }
3808: #else
3809: *file = NULL;
3810: *func = NULL;
3811: *line = 0;
3812: #endif
3813: PetscFunctionReturn(PETSC_SUCCESS);
3814: }
3816: /*@
3817: VecLockReadPush - Pushes a read-only lock on a vector to prevent it from being written to
3819: Logically Collective
3821: Input Parameter:
3822: . x - the vector
3824: Level: intermediate
3826: Notes:
3827: If this is set then calls to `VecGetArray()` or `VecSetValues()` or any other routines that change the vectors values will generate an error.
3829: The call can be nested, i.e., called multiple times on the same vector, but each `VecLockReadPush()` has to have one matching
3830: `VecLockReadPop()`, which removes the latest read-only lock.
3832: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPop()`, `VecLockGet()`
3833: @*/
3834: PetscErrorCode VecLockReadPush(Vec x)
3835: {
3836: PetscFunctionBegin;
3838: 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");
3839: #if !PetscDefined(HAVE_THREADSAFETY)
3840: {
3841: const char *file, *func;
3842: int index, line;
3844: if ((index = petscstack.currentsize - 2) == -1) {
3845: // vec was locked "outside" of petsc, either in user-land or main. the error message will
3846: // now show this function as the culprit, but it will include the stacktrace
3847: file = "unknown user-file";
3848: func = "unknown_user_function";
3849: line = 0;
3850: } else {
3851: PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected petscstack, have negative index %d", index);
3852: file = petscstack.file[index];
3853: func = petscstack.function[index];
3854: line = petscstack.line[index];
3855: }
3856: PetscStackPush_Private(x->lockstack, file, func, line, petscstack.petscroutine[index], PETSC_FALSE);
3857: }
3858: #endif
3859: PetscFunctionReturn(PETSC_SUCCESS);
3860: }
3862: /*@
3863: VecLockReadPop - Pops a read-only lock from a vector
3865: Logically Collective
3867: Input Parameter:
3868: . x - the vector
3870: Level: intermediate
3872: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockGet()`
3873: @*/
3874: PetscErrorCode VecLockReadPop(Vec x)
3875: {
3876: PetscFunctionBegin;
3878: PetscCheck(--x->lock >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Vector has been unlocked from read-only access too many times");
3879: #if !PetscDefined(HAVE_THREADSAFETY)
3880: {
3881: const char *previous = x->lockstack.function[x->lockstack.currentsize - 1];
3883: PetscStackPop_Private(x->lockstack, previous);
3884: }
3885: #endif
3886: PetscFunctionReturn(PETSC_SUCCESS);
3887: }
3889: /*@C
3890: VecLockWriteSet - Lock or unlock a vector for exclusive read/write access
3892: Logically Collective
3894: Input Parameters:
3895: + x - the vector
3896: - flg - `PETSC_TRUE` to lock the vector for exclusive read/write access; `PETSC_FALSE` to unlock it.
3898: Level: intermediate
3900: Notes:
3901: The function is useful in split-phase computations, which usually have a begin phase and an end phase.
3902: One can call `VecLockWriteSet`(x,`PETSC_TRUE`) in the begin phase to lock a vector for exclusive
3903: access, and call `VecLockWriteSet`(x,`PETSC_FALSE`) in the end phase to unlock the vector from exclusive
3904: access. In this way, one is ensured no other operations can access the vector in between. The code may like
3906: .vb
3907: VecGetArray(x,&xdata); // begin phase
3908: VecLockWriteSet(v,PETSC_TRUE);
3910: Other operations, which can not access x anymore (they can access xdata, of course)
3912: VecRestoreArray(x,&vdata); // end phase
3913: VecLockWriteSet(v,PETSC_FALSE);
3914: .ve
3916: The call can not be nested on the same vector, in other words, one can not call `VecLockWriteSet`(x,`PETSC_TRUE`)
3917: again before calling `VecLockWriteSet`(v,`PETSC_FALSE`).
3919: .seealso: [](ch_vectors), `Vec`, `VecRestoreArray()`, `VecGetArrayRead()`, `VecLockReadPush()`, `VecLockReadPop()`, `VecLockGet()`
3920: @*/
3921: PetscErrorCode VecLockWriteSet(Vec x, PetscBool flg)
3922: {
3923: PetscFunctionBegin;
3925: if (flg) {
3926: 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");
3927: 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");
3928: x->lock = -1;
3929: } else {
3930: 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");
3931: x->lock = 0;
3932: }
3933: PetscFunctionReturn(PETSC_SUCCESS);
3934: }
3936: /*@
3937: VecLockPush - Pushes a read-only lock on a vector to prevent it from being written to
3939: Level: deprecated
3941: .seealso: [](ch_vectors), `Vec`, `VecLockReadPush()`
3942: @*/
3943: PetscErrorCode VecLockPush(Vec x)
3944: {
3945: PetscFunctionBegin;
3946: PetscCall(VecLockReadPush(x));
3947: PetscFunctionReturn(PETSC_SUCCESS);
3948: }
3950: /*@
3951: VecLockPop - Pops a read-only lock from a vector
3953: Level: deprecated
3955: .seealso: [](ch_vectors), `Vec`, `VecLockReadPop()`
3956: @*/
3957: PetscErrorCode VecLockPop(Vec x)
3958: {
3959: PetscFunctionBegin;
3960: PetscCall(VecLockReadPop(x));
3961: PetscFunctionReturn(PETSC_SUCCESS);
3962: }
3964: #endif