Actual source code: vector.c
1: /*
2: Provides the interface functions for vector operations that do NOT have PetscScalar/PetscReal in the signature
3: These are the vector functions the user calls.
4: */
5: #include <petsc/private/vecimpl.h>
6: #include <petsc/private/deviceimpl.h>
8: /* Logging support */
9: PetscClassId VEC_CLASSID;
10: PetscLogEvent VEC_View, VEC_Max, VEC_Min, VEC_Dot, VEC_MDot, VEC_TDot;
11: PetscLogEvent VEC_Norm, VEC_Normalize, VEC_Scale, VEC_Shift, VEC_Copy, VEC_Set, VEC_AXPY, VEC_AYPX, VEC_WAXPY;
12: PetscLogEvent VEC_MTDot, VEC_MAXPY, VEC_Swap, VEC_AssemblyBegin, VEC_ScatterBegin, VEC_ScatterEnd;
13: PetscLogEvent VEC_AssemblyEnd, VEC_PointwiseMult, VEC_PointwiseDivide, VEC_Reciprocal, VEC_SetValues, VEC_Load, VEC_SetPreallocateCOO, VEC_SetValuesCOO;
14: PetscLogEvent VEC_SetRandom, VEC_ReduceArithmetic, VEC_ReduceCommunication, VEC_ReduceBegin, VEC_ReduceEnd, VEC_Ops;
15: PetscLogEvent VEC_DotNorm2, VEC_AXPBYPCZ;
16: PetscLogEvent VEC_ViennaCLCopyFromGPU, VEC_ViennaCLCopyToGPU;
17: PetscLogEvent VEC_CUDACopyFromGPU, VEC_CUDACopyToGPU;
18: PetscLogEvent VEC_HIPCopyFromGPU, VEC_HIPCopyToGPU;
20: /*@
21: VecStashGetInfo - Gets how many values are currently in the vector stash, i.e. need
22: to be communicated to other processors during the `VecAssemblyBegin()`/`VecAssemblyEnd()` process
24: Not Collective
26: Input Parameter:
27: . vec - the vector
29: Output Parameters:
30: + nstash - the size of the stash
31: . reallocs - the number of additional mallocs incurred in building the stash
32: . bnstash - the size of the block stash
33: - breallocs - the number of additional mallocs incurred in building the block stash (from `VecSetValuesBlocked()`)
35: Level: advanced
37: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecStashSetInitialSize()`, `VecStashView()`
38: @*/
39: PetscErrorCode VecStashGetInfo(Vec vec, PetscInt *nstash, PetscInt *reallocs, PetscInt *bnstash, PetscInt *breallocs)
40: {
41: PetscFunctionBegin;
42: PetscCall(VecStashGetInfo_Private(&vec->stash, nstash, reallocs));
43: PetscCall(VecStashGetInfo_Private(&vec->bstash, bnstash, breallocs));
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: /*@
48: VecSetLocalToGlobalMapping - Sets a local numbering to global numbering used
49: by the routine `VecSetValuesLocal()` to allow users to insert vector entries
50: using a local (per-processor) numbering.
52: Logically Collective
54: Input Parameters:
55: + x - vector
56: - mapping - mapping created with `ISLocalToGlobalMappingCreate()` or `ISLocalToGlobalMappingCreateIS()`
58: Level: intermediate
60: Notes:
61: All vectors obtained with `VecDuplicate()` from this vector inherit the same mapping.
63: Vectors obtained with `DMCreateGlobaVector()` will often have this attribute attached to the vector so this call is not needed
65: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesLocal()`,
66: `VecGetLocalToGlobalMapping()`, `VecSetValuesBlockedLocal()`
67: @*/
68: PetscErrorCode VecSetLocalToGlobalMapping(Vec x, ISLocalToGlobalMapping mapping)
69: {
70: PetscFunctionBegin;
73: if (x->ops->setlocaltoglobalmapping) PetscUseTypeMethod(x, setlocaltoglobalmapping, mapping);
74: else PetscCall(PetscLayoutSetISLocalToGlobalMapping(x->map, mapping));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: /*@
79: VecGetLocalToGlobalMapping - Gets the local-to-global numbering set by `VecSetLocalToGlobalMapping()`
81: Not Collective
83: Input Parameter:
84: . X - the vector
86: Output Parameter:
87: . mapping - the mapping
89: Level: advanced
91: .seealso: [](ch_vectors), `Vec`, `VecSetValuesLocal()`, `VecSetLocalToGlobalMapping()`
92: @*/
93: PetscErrorCode VecGetLocalToGlobalMapping(Vec X, ISLocalToGlobalMapping *mapping)
94: {
95: PetscFunctionBegin;
98: PetscAssertPointer(mapping, 2);
99: if (X->ops->getlocaltoglobalmapping) PetscUseTypeMethod(X, getlocaltoglobalmapping, mapping);
100: else *mapping = X->map->mapping;
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: /*@
105: VecAssemblyBegin - Begins assembling the vector; that is ensuring all the vector's entries are stored on the correct MPI process. This routine should
106: be called after completing all calls to `VecSetValues()`.
108: Collective
110: Input Parameter:
111: . vec - the vector
113: Level: beginner
115: .seealso: [](ch_vectors), `Vec`, `VecAssemblyEnd()`, `VecSetValues()`
116: @*/
117: PetscErrorCode VecAssemblyBegin(Vec vec)
118: {
119: PetscFunctionBegin;
122: PetscCall(VecStashViewFromOptions(vec, NULL, "-vec_view_stash"));
123: PetscCall(PetscLogEventBegin(VEC_AssemblyBegin, vec, 0, 0, 0));
124: PetscTryTypeMethod(vec, assemblybegin);
125: PetscCall(PetscLogEventEnd(VEC_AssemblyBegin, vec, 0, 0, 0));
126: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: /*@
131: VecAssemblyEnd - Completes assembling the vector. This routine should be called after `VecAssemblyBegin()`.
133: Collective
135: Input Parameter:
136: . vec - the vector
138: Options Database Keys:
139: + -vec_view [viewertype][:...] - Display the vector. See `VecViewFromOptions()`/`PetscObjectViewFromOptions()` for the possible arguments
140: - -vecstash_view [viewertype][:...] - Display the vector stash. See `VecStashViewFromOptions()`/`PetscObjectViewFromOptions()` for the possible arguments
142: Level: beginner
144: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecSetValues()`, `VecViewFromOptions()`, `VecStashViewFromOptions()`,
145: `PetscObjectViewFromOptions()`
146: @*/
147: PetscErrorCode VecAssemblyEnd(Vec vec)
148: {
149: PetscFunctionBegin;
151: PetscCall(PetscLogEventBegin(VEC_AssemblyEnd, vec, 0, 0, 0));
153: PetscTryTypeMethod(vec, assemblyend);
154: PetscCall(PetscLogEventEnd(VEC_AssemblyEnd, vec, 0, 0, 0));
155: PetscCall(VecViewFromOptions(vec, NULL, "-vec_view"));
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: /*@
160: VecSetPreallocationCOO - set preallocation for a vector using a coordinate format of the entries with global indices
162: Collective
164: Input Parameters:
165: + x - vector being preallocated
166: . ncoo - number of entries
167: - coo_i - entry indices
169: Level: beginner
171: Notes:
172: This and `VecSetValuesCOO()` provide an alternative API to using `VecSetValues()` to provide vector values.
174: This API is particularly efficient for use on GPUs.
176: Entries can be repeated, see `VecSetValuesCOO()`. Negative indices are not allowed unless vector option `VEC_IGNORE_NEGATIVE_INDICES` is set,
177: in which case they, along with the corresponding entries in `VecSetValuesCOO()`, are ignored. If vector option `VEC_NO_OFF_PROC_ENTRIES` is set,
178: remote entries are ignored, otherwise, they will be properly added or inserted to the vector.
180: The array coo_i[] may be freed immediately after calling this function.
182: .seealso: [](ch_vectors), `Vec`, `VecSetValuesCOO()`, `VecSetPreallocationCOOLocal()`
183: @*/
184: PetscErrorCode VecSetPreallocationCOO(Vec x, PetscCount ncoo, const PetscInt coo_i[])
185: {
186: PetscFunctionBegin;
189: if (ncoo) PetscAssertPointer(coo_i, 3);
190: PetscCall(PetscLogEventBegin(VEC_SetPreallocateCOO, x, 0, 0, 0));
191: PetscCall(PetscLayoutSetUp(x->map));
192: if (x->ops->setpreallocationcoo) {
193: PetscUseTypeMethod(x, setpreallocationcoo, ncoo, coo_i);
194: } else {
195: PetscInt ncoo_i;
196: IS is_coo_i;
198: PetscCall(PetscIntCast(ncoo, &ncoo_i));
199: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncoo_i, coo_i, PETSC_COPY_VALUES, &is_coo_i));
200: PetscCall(PetscObjectCompose((PetscObject)x, "__PETSc_coo_i", (PetscObject)is_coo_i));
201: PetscCall(ISDestroy(&is_coo_i));
202: }
203: PetscCall(PetscLogEventEnd(VEC_SetPreallocateCOO, x, 0, 0, 0));
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: /*@
208: VecSetPreallocationCOOLocal - set preallocation for vectors using a coordinate format of the entries with local indices
210: Collective
212: Input Parameters:
213: + x - vector being preallocated
214: . ncoo - number of entries
215: - coo_i - row indices (local numbering; may be modified)
217: Level: beginner
219: Notes:
220: This and `VecSetValuesCOO()` provide an alternative API to using `VecSetValuesLocal()` to provide vector values.
222: This API is particularly efficient for use on GPUs.
224: The local indices are translated using the local to global mapping, thus `VecSetLocalToGlobalMapping()` must have been
225: called prior to this function.
227: The indices coo_i may be modified within this function. They might be translated to corresponding global
228: indices, but the caller should not rely on them having any specific value after this function returns. The arrays
229: can be freed or reused immediately after this function returns.
231: Entries can be repeated. Negative indices and remote indices might be allowed. see `VecSetPreallocationCOO()`.
233: .seealso: [](ch_vectors), `Vec`, `VecSetPreallocationCOO()`, `VecSetValuesCOO()`
234: @*/
235: PetscErrorCode VecSetPreallocationCOOLocal(Vec x, PetscCount ncoo, PetscInt coo_i[])
236: {
237: PetscInt ncoo_i;
238: ISLocalToGlobalMapping ltog;
240: PetscFunctionBegin;
243: if (ncoo) PetscAssertPointer(coo_i, 3);
244: PetscCall(PetscIntCast(ncoo, &ncoo_i));
245: PetscCall(PetscLayoutSetUp(x->map));
246: PetscCall(VecGetLocalToGlobalMapping(x, <og));
247: if (ltog) PetscCall(ISLocalToGlobalMappingApply(ltog, ncoo_i, coo_i, coo_i));
248: PetscCall(VecSetPreallocationCOO(x, ncoo, coo_i));
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: /*@
253: VecSetValuesCOO - set values at once in a vector preallocated using `VecSetPreallocationCOO()`
255: Collective
257: Input Parameters:
258: + x - vector being set
259: . coo_v - the value array
260: - imode - the insert mode
262: Level: beginner
264: Note:
265: This and `VecSetPreallocationCOO() or ``VecSetPreallocationCOOLocal()` provide an alternative API to using `VecSetValues()` to provide vector values.
267: This API is particularly efficient for use on GPUs.
269: The values must follow the order of the indices prescribed with `VecSetPreallocationCOO()` or `VecSetPreallocationCOOLocal()`.
270: When repeated entries are specified in the COO indices the `coo_v` values are first properly summed, regardless of the value of `imode`.
271: The imode flag indicates if `coo_v` must be added to the current values of the vector (`ADD_VALUES`) or overwritten (`INSERT_VALUES`).
272: `VecAssemblyBegin()` and `VecAssemblyEnd()` do not need to be called after this routine. It automatically handles the assembly process.
274: .seealso: [](ch_vectors), `Vec`, `VecSetPreallocationCOO()`, `VecSetPreallocationCOOLocal()`, `VecSetValues()`
275: @*/
276: PetscErrorCode VecSetValuesCOO(Vec x, const PetscScalar coo_v[], InsertMode imode)
277: {
278: PetscFunctionBegin;
282: PetscCall(PetscLogEventBegin(VEC_SetValuesCOO, x, 0, 0, 0));
283: if (x->ops->setvaluescoo) {
284: PetscUseTypeMethod(x, setvaluescoo, coo_v, imode);
285: PetscCall(PetscObjectStateIncrease((PetscObject)x));
286: } else {
287: IS is_coo_i;
288: const PetscInt *coo_i;
289: PetscInt ncoo;
290: PetscMemType mtype;
292: PetscCall(PetscGetMemType(coo_v, &mtype));
293: PetscCheck(mtype == PETSC_MEMTYPE_HOST, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_WRONG, "The basic VecSetValuesCOO() only supports v[] on host");
294: PetscCall(PetscObjectQuery((PetscObject)x, "__PETSc_coo_i", (PetscObject *)&is_coo_i));
295: PetscCheck(is_coo_i, PetscObjectComm((PetscObject)x), PETSC_ERR_COR, "Missing coo_i IS");
296: PetscCall(ISGetLocalSize(is_coo_i, &ncoo));
297: PetscCall(ISGetIndices(is_coo_i, &coo_i));
298: if (imode != ADD_VALUES) PetscCall(VecZeroEntries(x));
299: PetscCall(VecSetValues(x, ncoo, coo_i, coo_v, ADD_VALUES));
300: PetscCall(ISRestoreIndices(is_coo_i, &coo_i));
301: PetscCall(VecAssemblyBegin(x));
302: PetscCall(VecAssemblyEnd(x));
303: }
304: PetscCall(PetscLogEventEnd(VEC_SetValuesCOO, x, 0, 0, 0));
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }
308: static PetscErrorCode VecPointwiseApply_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx, PetscLogEvent event, const char async_name[], PetscErrorCode (*const pointwise_op)(Vec, Vec, Vec))
309: {
310: PetscErrorCode (*async_fn)(Vec, Vec, Vec, PetscDeviceContext) = NULL;
312: PetscFunctionBegin;
319: PetscCheckSameTypeAndComm(x, 2, y, 3);
320: PetscCheckSameTypeAndComm(y, 3, w, 1);
321: VecCheckSameSize(w, 1, x, 2);
322: VecCheckSameSize(w, 1, y, 3);
323: VecCheckAssembled(x);
324: VecCheckAssembled(y);
325: PetscCall(VecSetErrorIfLocked(w, 1));
328: if (dctx) PetscCall(PetscObjectQueryFunction((PetscObject)w, async_name, &async_fn));
329: if (event) PetscCall(PetscLogEventBegin(event, x, y, w, 0));
330: if (async_fn) PetscCall((*async_fn)(w, x, y, dctx));
331: else PetscCall((*pointwise_op)(w, x, y));
332: if (event) PetscCall(PetscLogEventEnd(event, x, y, w, 0));
333: PetscCall(PetscObjectStateIncrease((PetscObject)w));
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: PetscErrorCode VecPointwiseMaxAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
338: {
339: PetscFunctionBegin;
340: // REVIEW ME: no log event?
341: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMax), w->ops->pointwisemax));
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: /*@
346: VecPointwiseMax - Computes the component-wise maximum `w[i] = max(x[i], y[i])`.
348: Logically Collective
350: Input Parameters:
351: + x - the first input vector
352: - y - the second input vector
354: Output Parameter:
355: . w - the result
357: Level: advanced
359: Notes:
360: Any subset of the `x`, `y`, and `w` may be the same vector.
362: For complex numbers compares only the real part
364: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
365: @*/
366: PetscErrorCode VecPointwiseMax(Vec w, Vec x, Vec y)
367: {
368: PetscFunctionBegin;
369: PetscCall(VecPointwiseMaxAsync_Private(w, x, y, NULL));
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: PetscErrorCode VecPointwiseMinAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
374: {
375: PetscFunctionBegin;
376: // REVIEW ME: no log event?
377: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMin), w->ops->pointwisemin));
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: /*@
382: VecPointwiseMin - Computes the component-wise minimum `w[i] = min(x[i], y[i])`.
384: Logically Collective
386: Input Parameters:
387: + x - the first input vector
388: - y - the second input vector
390: Output Parameter:
391: . w - the result
393: Level: advanced
395: Notes:
396: Any subset of the `x`, `y`, and `w` may be the same vector.
398: For complex numbers compares only the real part
400: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
401: @*/
402: PetscErrorCode VecPointwiseMin(Vec w, Vec x, Vec y)
403: {
404: PetscFunctionBegin;
405: PetscCall(VecPointwiseMinAsync_Private(w, x, y, NULL));
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: PetscErrorCode VecPointwiseMaxAbsAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
410: {
411: PetscFunctionBegin;
412: // REVIEW ME: no log event?
413: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMaxAbs), w->ops->pointwisemaxabs));
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: /*@
418: VecPointwiseMaxAbs - Computes the component-wise maximum of the absolute values `w[i] = max(abs(x[i]), abs(y[i]))`.
420: Logically Collective
422: Input Parameters:
423: + x - the first input vector
424: - y - the second input vector
426: Output Parameter:
427: . w - the result
429: Level: advanced
431: Notes:
432: Any subset of the `x`, `y`, and `w` may be the same vector.
434: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMin()`, `VecPointwiseMax()`, `VecMaxPointwiseDivide()`
435: @*/
436: PetscErrorCode VecPointwiseMaxAbs(Vec w, Vec x, Vec y)
437: {
438: PetscFunctionBegin;
439: PetscCall(VecPointwiseMaxAbsAsync_Private(w, x, y, NULL));
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: PetscErrorCode VecPointwiseDivideAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
444: {
445: PetscFunctionBegin;
446: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, VEC_PointwiseDivide, VecAsyncFnName(PointwiseDivide), w->ops->pointwisedivide));
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: /*@
451: VecPointwiseDivide - Computes the component-wise division `w[i] = x[i] / y[i]`.
453: Logically Collective
455: Input Parameters:
456: + x - the numerator vector
457: - y - the denominator vector
459: Output Parameter:
460: . w - the result
462: Level: advanced
464: Note:
465: Any subset of the `x`, `y`, and `w` may be the same vector.
467: .seealso: [](ch_vectors), `Vec`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
468: @*/
469: PetscErrorCode VecPointwiseDivide(Vec w, Vec x, Vec y)
470: {
471: PetscFunctionBegin;
472: PetscCall(VecPointwiseDivideAsync_Private(w, x, y, NULL));
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: #define VEC_POINTWISE_SIGN_LOOP(y, x, n, func) \
477: PetscPragmaSIMD \
478: for (PetscInt i = 0; i < (n); i++) (y)[i] = func(PetscRealPart((x)[i]))
480: #define VEC_POINTWISE_SIGN_DISPATCH(y, x, n, sign_type) \
481: do { \
482: switch (sign_type) { \
483: case VEC_SIGN_ZERO_TO_ZERO: \
484: VEC_POINTWISE_SIGN_LOOP(y, x, n, VecSignZeroToZero_Private); \
485: break; \
486: case VEC_SIGN_ZERO_TO_SIGNED_ZERO: \
487: VEC_POINTWISE_SIGN_LOOP(y, x, n, VecSignZeroToSignedZero_Private); \
488: break; \
489: case VEC_SIGN_ZERO_TO_SIGNED_UNIT: \
490: VEC_POINTWISE_SIGN_LOOP(y, x, n, VecSignZeroToSignedUnit_Private); \
491: break; \
492: default: \
493: PetscUnreachable(); \
494: } \
495: } while (0)
497: PetscErrorCode VecPointwiseSignAsync_Private(Vec y, Vec x, VecSignMode sign_type, PetscDeviceContext dctx)
498: {
499: PetscOffloadMask mask;
500: PetscBool is_host;
501: PetscErrorCode (*async_fn)(Vec, Vec, VecSignMode, PetscDeviceContext) = NULL;
503: PetscFunctionBegin;
508: VecCheckSameSize(y, 1, x, 2);
509: VecCheckAssembled(x);
510: VecCheckAssembled(y);
511: PetscCall(VecSetErrorIfLocked(y, 1));
513: PetscCall(VecGetOffloadMask(x, &mask));
514: is_host = PetscOffloadHost(mask) ? PETSC_TRUE : PETSC_FALSE;
515: if (!is_host) PetscCall(PetscObjectQueryFunction((PetscObject)y, VEC_ASYNC_FN_NAME("PointwiseSign"), &async_fn));
516: if (async_fn) PetscCall((*async_fn)(y, x, sign_type, dctx));
517: else {
518: PetscInt n;
520: PetscCall(VecGetLocalSize(y, &n));
521: if (y == x) {
522: PetscScalar *_y;
524: PetscCall(VecGetArray(y, &_y));
525: VEC_POINTWISE_SIGN_DISPATCH(_y, _y, n, sign_type);
526: PetscCall(VecRestoreArray(y, &_y));
527: } else {
528: PetscScalar *_y;
529: const PetscScalar *_x;
531: PetscCall(VecGetArrayWrite(y, &_y));
532: PetscCall(VecGetArrayRead(x, &_x));
533: VEC_POINTWISE_SIGN_DISPATCH(_y, _x, n, sign_type);
534: PetscCall(VecRestoreArrayRead(x, &_x));
535: PetscCall(VecRestoreArrayWrite(y, &_y));
536: }
537: }
538: PetscCall(PetscObjectStateIncrease((PetscObject)y));
539: PetscFunctionReturn(PETSC_SUCCESS);
540: }
542: /*@
543: VecPointwiseSign - Computes the component-wise sign `y[i] = sign(x[i])`.
545: Logically Collective
547: Input Parameters:
548: + x - the input vector
549: - sign_type - `VecSignMode` indicating how the function should map zero values.
551: Output Parameter:
552: . y - the sign vector of `x`
554: Level: beginner
556: .seealso: [](ch_vectors), `Vec`, `VecSignMode`
557: @*/
558: PetscErrorCode VecPointwiseSign(Vec y, Vec x, VecSignMode sign_type)
559: {
560: PetscFunctionBegin;
561: PetscCall(VecPointwiseSignAsync_Private(y, x, sign_type, NULL));
562: PetscFunctionReturn(PETSC_SUCCESS);
563: }
565: PetscErrorCode VecPointwiseMultAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
566: {
567: PetscFunctionBegin;
569: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, VEC_PointwiseMult, VecAsyncFnName(PointwiseMult), w->ops->pointwisemult));
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: /*@
574: VecPointwiseMult - Computes the component-wise multiplication `w[i] = x[i] * y[i]`.
576: Logically Collective
578: Input Parameters:
579: + x - the first vector
580: - y - the second vector
582: Output Parameter:
583: . w - the result
585: Level: advanced
587: Note:
588: Any subset of the `x`, `y`, and `w` may be the same vector.
590: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
591: @*/
592: PetscErrorCode VecPointwiseMult(Vec w, Vec x, Vec y)
593: {
594: PetscFunctionBegin;
595: PetscCall(VecPointwiseMultAsync_Private(w, x, y, NULL));
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
599: /*@
600: VecDuplicate - Creates a new vector of the same type as an existing vector.
602: Collective
604: Input Parameter:
605: . v - a vector to mimic
607: Output Parameter:
608: . newv - location to put new vector
610: Level: beginner
612: Notes:
613: `VecDuplicate()` DOES NOT COPY the vector entries, but rather allocates storage
614: for the new vector. Use `VecCopy()` to copy a vector.
616: PETSc `Vec` always have all zero entries when created with `VecDuplicate()` until routines such as `VecSet()` or `VecSetValues()`
617: are used to change the values. There is no reason to call `VecZeroEntries()` after creation.
619: Use `VecDestroy()` to free the space. Use `VecDuplicateVecs()` to get several
620: vectors.
622: .seealso: [](ch_vectors), `Vec`, `VecDestroy()`, `VecDuplicateVecs()`, `VecCreate()`, `VecCopy()`
623: @*/
624: PetscErrorCode VecDuplicate(Vec v, Vec *newv)
625: {
626: PetscFunctionBegin;
628: PetscAssertPointer(newv, 2);
630: PetscUseTypeMethod(v, duplicate, newv);
631: #if PetscDefined(HAVE_DEVICE)
632: if (v->boundtocpu && v->bindingpropagates) {
633: PetscCall(VecSetBindingPropagates(*newv, PETSC_TRUE));
634: PetscCall(VecBindToCPU(*newv, PETSC_TRUE));
635: }
636: #endif
637: PetscCall(PetscObjectStateIncrease((PetscObject)*newv));
638: PetscFunctionReturn(PETSC_SUCCESS);
639: }
641: /*@
642: VecDestroy - Destroys a vector.
644: Collective
646: Input Parameter:
647: . v - the vector
649: Level: beginner
651: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecDuplicate()`, `VecDestroyVecs()`
652: @*/
653: PetscErrorCode VecDestroy(Vec *v)
654: {
655: PetscFunctionBegin;
656: PetscAssertPointer(v, 1);
657: if (!*v) PetscFunctionReturn(PETSC_SUCCESS);
659: if (--((PetscObject)*v)->refct > 0) {
660: *v = NULL;
661: PetscFunctionReturn(PETSC_SUCCESS);
662: }
664: PetscCall(PetscObjectSAWsViewOff((PetscObject)*v));
665: /* destroy the internal part */
666: PetscTryTypeMethod(*v, destroy);
667: PetscCall(PetscFree((*v)->defaultrandtype));
668: /* destroy the external/common part */
669: PetscCall(PetscLayoutDestroy(&(*v)->map));
670: PetscCall(PetscHeaderDestroy(v));
671: PetscFunctionReturn(PETSC_SUCCESS);
672: }
674: /*@C
675: VecDuplicateVecs - Creates several vectors of the same type as an existing vector.
677: Collective
679: Input Parameters:
680: + m - the number of vectors to obtain
681: - v - a vector to mimic
683: Output Parameter:
684: . V - location to put pointer to array of vectors
686: Level: intermediate
688: Notes:
689: Use `VecDestroyVecs()` to free the space. Use `VecDuplicate()` to form a single
690: vector.
692: PETSc `Vec` always have all zero entries when created with `VecDuplicateVecs()` until routines such as `VecSet()` or `VecSetValues()`
693: are used to change the values. There is no reason to call `VecZeroEntries()` after creation.
695: Some implementations ensure that the arrays accessed by each vector are contiguous in memory. Certain `VecMDot()` and `VecMAXPY()`
696: implementations utilize this property to use BLAS 2 operations for higher efficiency. This is especially useful in `KSPGMRES`, see
697: `KSPGMRESSetPreAllocateVectors()`.
699: Fortran Note:
700: .vb
701: Vec, pointer :: V(:)
702: .ve
704: .seealso: [](ch_vectors), `Vec`, [](ch_fortran), `VecDestroyVecs()`, `VecDuplicate()`, `VecCreate()`, `VecMDot()`, `VecMAXPY()`, `KSPGMRES`,
705: `KSPGMRESSetPreAllocateVectors()`
706: @*/
707: PetscErrorCode VecDuplicateVecs(Vec v, PetscInt m, Vec *V[])
708: {
709: PetscFunctionBegin;
711: PetscAssertPointer(V, 3);
713: PetscUseTypeMethod(v, duplicatevecs, m, V);
714: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
715: if (v->boundtocpu && v->bindingpropagates) {
716: PetscInt i;
718: for (i = 0; i < m; i++) {
719: /* Since ops->duplicatevecs might itself propagate the value of boundtocpu,
720: * avoid unnecessary overhead by only calling VecBindToCPU() if the vector isn't already bound. */
721: if (!(*V)[i]->boundtocpu) {
722: PetscCall(VecSetBindingPropagates((*V)[i], PETSC_TRUE));
723: PetscCall(VecBindToCPU((*V)[i], PETSC_TRUE));
724: }
725: }
726: }
727: #endif
728: PetscFunctionReturn(PETSC_SUCCESS);
729: }
731: /*@C
732: VecDestroyVecs - Frees a block of vectors obtained with `VecDuplicateVecs()`.
734: Collective
736: Input Parameters:
737: + m - the number of vectors previously obtained, if zero no vectors are destroyed
738: - vv - pointer to pointer to array of vector pointers, if `NULL` no vectors are destroyed
740: Level: intermediate
742: .seealso: [](ch_vectors), `Vec`, [](ch_fortran), `VecDuplicateVecs()`, `VecDestroyVecsf90()`
743: @*/
744: PetscErrorCode VecDestroyVecs(PetscInt m, Vec *vv[])
745: {
746: PetscFunctionBegin;
747: PetscAssertPointer(vv, 2);
748: PetscCheck(m >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Trying to destroy negative number of vectors %" PetscInt_FMT, m);
749: if (!m || !*vv) {
750: *vv = NULL;
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
755: PetscCall((*(**vv)->ops->destroyvecs)(m, *vv));
756: *vv = NULL;
757: PetscFunctionReturn(PETSC_SUCCESS);
758: }
760: /*@
761: VecViewFromOptions - View a vector based on values in the options database
763: Collective
765: Input Parameters:
766: + A - the vector
767: . obj - optional object that provides the options prefix for this viewing, use 'NULL' to use the prefix of `A`
768: - name - command line option
770: Options Database Key:
771: . -name [viewertype][:...] - option name and values. See `PetscObjectViewFromOptions()` for the possible arguments
773: Level: intermediate
775: .seealso: [](ch_vectors), `Vec`, `VecView`, `PetscObjectViewFromOptions()`, `VecCreate()`
776: @*/
777: PetscErrorCode VecViewFromOptions(Vec A, PeOp PetscObject obj, const char name[])
778: {
779: PetscFunctionBegin;
781: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
782: PetscFunctionReturn(PETSC_SUCCESS);
783: }
785: /*@
786: VecView - Views a vector object.
788: Collective
790: Input Parameters:
791: + vec - the vector
792: - viewer - an optional `PetscViewer` visualization context
794: Level: beginner
796: Notes:
797: The available visualization contexts include
798: + `PETSC_VIEWER_STDOUT_SELF` - for sequential vectors
799: . `PETSC_VIEWER_STDOUT_WORLD` - for parallel vectors created on `PETSC_COMM_WORLD`
800: - `PETSC_VIEWER_STDOUT`_(comm) - for parallel vectors created on MPI communicator comm
802: You can change the format the vector is printed using the
803: option `PetscViewerPushFormat()`.
805: The user can open alternative viewers with
806: + `PetscViewerASCIIOpen()` - Outputs vector to a specified file
807: . `PetscViewerBinaryOpen()` - Outputs vector in binary to a
808: specified file; corresponding input uses `VecLoad()`
809: . `PetscViewerDrawOpen()` - Outputs vector to an X window display
810: . `PetscViewerSocketOpen()` - Outputs vector to Socket viewer
811: - `PetscViewerHDF5Open()` - Outputs vector to HDF5 file viewer
813: The user can call `PetscViewerPushFormat()` to specify the output
814: format of ASCII printed objects (when using `PETSC_VIEWER_STDOUT_SELF`,
815: `PETSC_VIEWER_STDOUT_WORLD` and `PetscViewerASCIIOpen()`). Available formats include
816: + `PETSC_VIEWER_DEFAULT` - default, prints vector contents
817: . `PETSC_VIEWER_ASCII_MATLAB` - prints vector contents in MATLAB format
818: . `PETSC_VIEWER_ASCII_INDEX` - prints vector contents, including indices of vector elements
819: - `PETSC_VIEWER_ASCII_COMMON` - prints vector contents, using a
820: format common among all vector types
822: You can pass any number of vector objects, or other PETSc objects to the same viewer.
824: In the debugger you can do call `VecView`(v,0) to display the vector. (The same holds for any PETSc object viewer).
826: Notes for binary viewer:
827: If you pass multiple vectors to a binary viewer you can read them back in the same order
828: with `VecLoad()`.
830: If the blocksize of the vector is greater than one then you must provide a unique prefix to
831: the vector with `PetscObjectSetOptionsPrefix`((`PetscObject`)vec,"uniqueprefix"); BEFORE calling `VecView()` on the
832: vector to be stored and then set that same unique prefix on the vector that you pass to `VecLoad()`. The blocksize
833: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
834: filename. If you copy the binary file, make sure you copy the associated .info file with it.
836: See the manual page for `VecLoad()` on the exact format the binary viewer stores
837: the values in the file.
839: Notes for HDF5 Viewer:
840: The name of the `Vec` (given with `PetscObjectSetName()` is the name that is used
841: for the object in the HDF5 file. If you wish to store the same Vec into multiple
842: datasets in the same file (typically with different values), you must change its
843: name each time before calling the `VecView()`. To load the same vector,
844: the name of the Vec object passed to `VecLoad()` must be the same.
846: If the block size of the vector is greater than 1 then it is used as the first dimension in the HDF5 array.
847: If the function `PetscViewerHDF5SetBaseDimension2()`is called then even if the block size is one it will
848: be used as the first dimension in the HDF5 array (that is the HDF5 array will always be two dimensional)
849: See also `PetscViewerHDF5SetTimestep()` which adds an additional complication to reading and writing `Vec`
850: with the HDF5 viewer.
852: .seealso: [](ch_vectors), `Vec`, `VecViewFromOptions()`, `PetscViewerASCIIOpen()`, `PetscViewerDrawOpen()`, `PetscDrawLGCreate()`,
853: `PetscViewerSocketOpen()`, `PetscViewerBinaryOpen()`, `VecLoad()`, `PetscViewerCreate()`,
854: `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`, `PetscViewerHDF5SetTimestep()`
855: @*/
856: PetscErrorCode VecView(Vec vec, PetscViewer viewer)
857: {
858: PetscBool isascii;
859: PetscViewerFormat format;
860: PetscMPIInt size;
862: PetscFunctionBegin;
865: VecCheckAssembled(vec);
866: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec), &viewer));
868: PetscCall(PetscViewerGetFormat(viewer, &format));
869: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
870: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) PetscFunctionReturn(PETSC_SUCCESS);
872: PetscCheck(!vec->stash.n && !vec->bstash.n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call VecAssemblyBegin/End() before viewing this vector");
874: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
875: if (isascii) {
876: PetscInt rows, bs;
878: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)vec, viewer));
879: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
880: PetscCall(PetscViewerASCIIPushTab(viewer));
881: PetscCall(VecGetSize(vec, &rows));
882: PetscCall(VecGetBlockSize(vec, &bs));
883: if (bs != 1) {
884: PetscCall(PetscViewerASCIIPrintf(viewer, "length=%" PetscInt_FMT ", bs=%" PetscInt_FMT "\n", rows, bs));
885: } else {
886: PetscCall(PetscViewerASCIIPrintf(viewer, "length=%" PetscInt_FMT "\n", rows));
887: }
888: PetscCall(PetscViewerASCIIPopTab(viewer));
889: }
890: }
891: PetscCall(VecLockReadPush(vec));
892: PetscCall(PetscLogEventBegin(VEC_View, vec, viewer, 0, 0));
893: if ((format == PETSC_VIEWER_NATIVE || format == PETSC_VIEWER_LOAD_BALANCE) && vec->ops->viewnative) {
894: PetscUseTypeMethod(vec, viewnative, viewer);
895: } else {
896: PetscUseTypeMethod(vec, view, viewer);
897: }
898: PetscCall(VecLockReadPop(vec));
899: PetscCall(PetscLogEventEnd(VEC_View, vec, viewer, 0, 0));
900: PetscFunctionReturn(PETSC_SUCCESS);
901: }
903: #if defined(PETSC_USE_DEBUG)
904: #include <../src/sys/totalview/tv_data_display.h>
905: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
906: {
907: const PetscScalar *values;
908: char type[32];
910: TV_add_row("Local rows", "int", &v->map->n);
911: TV_add_row("Global rows", "int", &v->map->N);
912: TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
913: PetscCall(VecGetArrayRead((Vec)v, &values));
914: PetscCall(PetscSNPrintf(type, 32, "double[%" PetscInt_FMT "]", v->map->n));
915: TV_add_row("values", type, values);
916: PetscCall(VecRestoreArrayRead((Vec)v, &values));
917: return TV_format_OK;
918: }
919: #endif
921: /*@C
922: VecViewNative - Views a vector object with the original type specific viewer
924: Collective
926: Input Parameters:
927: + vec - the vector
928: - viewer - an optional `PetscViewer` visualization context
930: Level: developer
932: Note:
933: This can be used with, for example, vectors obtained with `DMCreateGlobalVector()` for a `DMDA` to display the vector
934: in the PETSc storage format (each MPI process values follow the previous MPI processes) instead of the "natural" grid
935: ordering.
937: .seealso: [](ch_vectors), `Vec`, `PetscViewerASCIIOpen()`, `PetscViewerDrawOpen()`, `PetscDrawLGCreate()`, `VecView()`,
938: `PetscViewerSocketOpen()`, `PetscViewerBinaryOpen()`, `VecLoad()`, `PetscViewerCreate()`,
939: `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`, `PetscViewerHDF5SetTimestep()`
940: @*/
941: PetscErrorCode VecViewNative(Vec vec, PetscViewer viewer)
942: {
943: PetscFunctionBegin;
946: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec), &viewer));
948: PetscUseTypeMethod(vec, viewnative, viewer);
949: PetscFunctionReturn(PETSC_SUCCESS);
950: }
952: /*@
953: VecGetSize - Returns the global number of elements of the vector.
955: Not Collective
957: Input Parameter:
958: . x - the vector
960: Output Parameter:
961: . size - the global length of the vector
963: Level: beginner
965: .seealso: [](ch_vectors), `Vec`, `VecGetLocalSize()`
966: @*/
967: PetscErrorCode VecGetSize(Vec x, PetscInt *size)
968: {
969: PetscFunctionBegin;
971: PetscAssertPointer(size, 2);
973: PetscUseTypeMethod(x, getsize, size);
974: PetscFunctionReturn(PETSC_SUCCESS);
975: }
977: /*@
978: VecGetLocalSize - Returns the number of elements of the vector stored
979: in local memory (that is on this MPI process)
981: Not Collective
983: Input Parameter:
984: . x - the vector
986: Output Parameter:
987: . size - the length of the local piece of the vector
989: Level: beginner
991: .seealso: [](ch_vectors), `Vec`, `VecGetSize()`
992: @*/
993: PetscErrorCode VecGetLocalSize(Vec x, PetscInt *size)
994: {
995: PetscFunctionBegin;
997: PetscAssertPointer(size, 2);
999: PetscUseTypeMethod(x, getlocalsize, size);
1000: PetscFunctionReturn(PETSC_SUCCESS);
1001: }
1003: /*@
1004: VecGetOwnershipRange - Returns the range of indices owned by
1005: this process. The vector is laid out with the
1006: first `n1` elements on the first processor, next `n2` elements on the
1007: second, etc. For certain parallel layouts this range may not be
1008: well defined.
1010: Not Collective
1012: Input Parameter:
1013: . x - the vector
1015: Output Parameters:
1016: + low - the first local element, pass in `NULL` if not interested
1017: - high - one more than the last local element, pass in `NULL` if not interested
1019: Level: beginner
1021: Notes:
1022: If the `Vec` was obtained from a `DM` with `DMCreateGlobalVector()`, then the range values are determined by the specific `DM`.
1024: If the `Vec` was created directly the range values are determined by the local size passed to `VecSetSizes()` or `VecCreateMPI()`.
1025: If `PETSC_DECIDE` was passed as the local size, then the vector uses default values for the range using `PetscSplitOwnership()`.
1027: The high argument is one more than the last element stored locally.
1029: For certain `DM`, such as `DMDA`, it is better to use `DM` specific routines, such as `DMDAGetGhostCorners()`, to determine
1030: the local values in the vector.
1032: .seealso: [](ch_vectors), `Vec`, `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `VecGetOwnershipRanges()`, `PetscSplitOwnership()`,
1033: `VecSetSizes()`, `VecCreateMPI()`, `PetscLayout`, `DMDAGetGhostCorners()`, `DM`
1034: @*/
1035: PetscErrorCode VecGetOwnershipRange(Vec x, PetscInt *low, PetscInt *high)
1036: {
1037: PetscFunctionBegin;
1040: if (low) PetscAssertPointer(low, 2);
1041: if (high) PetscAssertPointer(high, 3);
1042: if (low) *low = x->map->rstart;
1043: if (high) *high = x->map->rend;
1044: PetscFunctionReturn(PETSC_SUCCESS);
1045: }
1047: /*@C
1048: VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
1049: The vector is laid out with the
1050: first `n1` elements on the first processor, next `n2` elements on the
1051: second, etc. For certain parallel layouts this range may not be
1052: well defined.
1054: Not Collective
1056: Input Parameter:
1057: . x - the vector
1059: Output Parameter:
1060: . ranges - array of length `size` + 1 with the start and end+1 for each process
1062: Level: beginner
1064: Notes:
1065: If the `Vec` was obtained from a `DM` with `DMCreateGlobalVector()`, then the range values are determined by the specific `DM`.
1067: If the `Vec` was created directly the range values are determined by the local size passed to `VecSetSizes()` or `VecCreateMPI()`.
1068: If `PETSC_DECIDE` was passed as the local size, then the vector uses default values for the range using `PetscSplitOwnership()`.
1070: The high argument is one more than the last element stored locally.
1072: For certain `DM`, such as `DMDA`, it is better to use `DM` specific routines, such as `DMDAGetGhostCorners()`, to determine
1073: the local values in the vector.
1075: The high argument is one more than the last element stored locally.
1077: If `ranges` are used after all vectors that share the ranges has been destroyed, then the program will crash accessing `ranges`.
1079: Fortran Note:
1080: The argument `ranges` must be declared as
1081: .vb
1082: PetscInt, pointer :: ranges(:)
1083: .ve
1084: and you have to return it with a call to `VecRestoreOwnershipRanges()` when no longer needed
1086: .seealso: [](ch_vectors), `Vec`, `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `VecGetOwnershipRange()`, `PetscSplitOwnership()`,
1087: `VecSetSizes()`, `VecCreateMPI()`, `PetscLayout`, `DMDAGetGhostCorners()`, `DM`
1088: @*/
1089: PetscErrorCode VecGetOwnershipRanges(Vec x, const PetscInt *ranges[])
1090: {
1091: PetscFunctionBegin;
1094: PetscCall(PetscLayoutGetRanges(x->map, ranges));
1095: PetscFunctionReturn(PETSC_SUCCESS);
1096: }
1098: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
1099: /*@
1100: VecSetOption - Sets an option for controlling a vector's behavior.
1102: Collective
1104: Input Parameters:
1105: + x - the vector
1106: . op - the option
1107: - flag - turn the option on or off
1109: Supported Options:
1110: + `VEC_IGNORE_OFF_PROC_ENTRIES` - which causes `VecSetValues()` to ignore
1111: entries destined to be stored on a separate processor. This can be used
1112: to eliminate the global reduction in the `VecAssemblyBegin()` if you know
1113: that you have only used `VecSetValues()` to set local elements
1114: . `VEC_IGNORE_NEGATIVE_INDICES` - which means you can pass negative indices
1115: in ix in calls to `VecSetValues()` or `VecGetValues()`. These rows are simply
1116: ignored.
1117: - `VEC_SUBSET_OFF_PROC_ENTRIES` - which causes `VecAssemblyBegin()` to assume that the off-process
1118: entries will always be a subset (possibly equal) of the off-process entries set on the
1119: first assembly which had a true `VEC_SUBSET_OFF_PROC_ENTRIES` and the vector has not
1120: changed this flag afterwards. If this assembly is not such first assembly, then this
1121: assembly can reuse the communication pattern setup in that first assembly, thus avoiding
1122: a global reduction. Subsequent assemblies setting off-process values should use the same
1123: InsertMode as the first assembly.
1125: Level: intermediate
1127: Developer Notes:
1128: The `InsertMode` restriction could be removed by packing the stash messages out of place.
1130: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`
1131: @*/
1132: PetscErrorCode VecSetOption(Vec x, VecOption op, PetscBool flag)
1133: {
1134: PetscFunctionBegin;
1137: PetscTryTypeMethod(x, setoption, op, flag);
1138: PetscFunctionReturn(PETSC_SUCCESS);
1139: }
1141: /* Default routines for obtaining and releasing; */
1142: /* may be used by any implementation */
1143: PetscErrorCode VecDuplicateVecs_Default(Vec w, PetscInt m, Vec *V[])
1144: {
1145: PetscFunctionBegin;
1146: PetscCheck(m > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "m must be > 0: m = %" PetscInt_FMT, m);
1147: PetscCall(PetscMalloc1(m, V));
1148: for (PetscInt i = 0; i < m; i++) PetscCall(VecDuplicate(w, *V + i));
1149: PetscFunctionReturn(PETSC_SUCCESS);
1150: }
1152: PetscErrorCode VecDestroyVecs_Default(PetscInt m, Vec v[])
1153: {
1154: PetscInt i;
1156: PetscFunctionBegin;
1157: PetscAssertPointer(v, 2);
1158: for (i = 0; i < m; i++) PetscCall(VecDestroy(&v[i]));
1159: PetscCall(PetscFree(v));
1160: PetscFunctionReturn(PETSC_SUCCESS);
1161: }
1163: /*@
1164: VecResetArray - Resets a vector to use its default memory. Call this
1165: after the use of `VecPlaceArray()`.
1167: Not Collective
1169: Input Parameter:
1170: . vec - the vector
1172: Level: developer
1174: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecPlaceArray()`
1175: @*/
1176: PetscErrorCode VecResetArray(Vec vec)
1177: {
1178: PetscFunctionBegin;
1181: PetscUseTypeMethod(vec, resetarray);
1182: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
1183: PetscFunctionReturn(PETSC_SUCCESS);
1184: }
1186: /*@
1187: VecLoad - Loads a vector that has been stored in binary or HDF5 format
1188: with `VecView()`.
1190: Collective
1192: Input Parameters:
1193: + vec - the newly loaded vector, this needs to have been created with `VecCreate()` or
1194: some related function before the call to `VecLoad()`.
1195: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()` or
1196: HDF5 file viewer, obtained from `PetscViewerHDF5Open()`
1198: Level: intermediate
1200: Notes:
1201: Defaults to the standard `VECSEQ` or `VECMPI`, if you want some other type of `Vec` call `VecSetFromOptions()`
1202: before calling this.
1204: The input file must contain the full global vector, as
1205: written by the routine `VecView()`.
1207: If the type or size of `vec` is not set before a call to `VecLoad()`, PETSc
1208: sets the type and the local and global sizes based on the vector it is reading in. If type and/or
1209: sizes are already set, then the same are used.
1211: If using the binary viewer and the blocksize of the vector is greater than one then you must provide a unique prefix to
1212: the vector with `PetscObjectSetOptionsPrefix`((`PetscObject`)vec,"uniqueprefix"); BEFORE calling `VecView()` on the
1213: vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
1214: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
1215: filename. If you copy the binary file, make sure you copy the associated .info file with it.
1217: If using HDF5, you must assign the `Vec` the same name as was used in the Vec
1218: that was stored in the file using `PetscObjectSetName()`. Otherwise you will
1219: get the error message: "Cannot H5DOpen2() with `Vec` name NAMEOFOBJECT".
1221: If the HDF5 file contains a two dimensional array the first dimension is treated as the block size
1222: in loading the vector. Hence, for example, using MATLAB notation h5create('vector.dat','/Test_Vec',[27 1]);
1223: will load a vector of size 27 and block size 27 thus resulting in all 27 entries being on the first process of
1224: vectors communicator and the rest of the processes having zero entries
1226: Notes for advanced users when using the binary viewer:
1227: Most users should not need to know the details of the binary storage
1228: format, since `VecLoad()` and `VecView()` completely hide these details.
1229: But for anyone who's interested, the standard binary vector storage
1230: format is
1231: .vb
1232: PetscInt VEC_FILE_CLASSID
1233: PetscInt number of rows
1234: PetscScalar *values of all entries
1235: .ve
1237: In addition, PETSc automatically uses byte swapping to work on all machines; the files
1238: are written ALWAYS using big-endian ordering. On small-endian machines the numbers
1239: are converted to the small-endian format when they are read in from the file.
1240: See PetscBinaryRead() and PetscBinaryWrite() to see how this may be done.
1242: .seealso: [](ch_vectors), `Vec`, `PetscViewerBinaryOpen()`, `VecView()`, `MatLoad()`
1243: @*/
1244: PetscErrorCode VecLoad(Vec vec, PetscViewer viewer)
1245: {
1246: PetscViewerFormat format;
1248: PetscFunctionBegin;
1251: PetscCheckSameComm(vec, 1, viewer, 2);
1253: PetscCall(VecSetErrorIfLocked(vec, 1));
1254: if (!((PetscObject)vec)->type_name && !vec->ops->create) PetscCall(VecSetType(vec, VECSTANDARD));
1255: PetscCall(PetscLogEventBegin(VEC_Load, viewer, 0, 0, 0));
1256: PetscCall(PetscViewerGetFormat(viewer, &format));
1257: if (format == PETSC_VIEWER_NATIVE && vec->ops->loadnative) {
1258: PetscUseTypeMethod(vec, loadnative, viewer);
1259: } else {
1260: PetscUseTypeMethod(vec, load, viewer);
1261: }
1262: PetscCall(PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0));
1263: PetscFunctionReturn(PETSC_SUCCESS);
1264: }
1266: /*@
1267: VecReciprocal - Replaces each component of a vector by its reciprocal.
1269: Logically Collective
1271: Input Parameter:
1272: . vec - the vector
1274: Output Parameter:
1275: . vec - the vector reciprocal
1277: Level: intermediate
1279: Note:
1280: Vector entries with value 0.0 are not changed
1282: .seealso: [](ch_vectors), `Vec`, `VecLog()`, `VecExp()`, `VecSqrtAbs()`
1283: @*/
1284: PetscErrorCode VecReciprocal(Vec vec)
1285: {
1286: PetscFunctionBegin;
1287: PetscCall(VecReciprocalAsync_Private(vec, NULL));
1288: PetscFunctionReturn(PETSC_SUCCESS);
1289: }
1291: /*@C
1292: VecSetOperation - Allows the user to override a particular vector operation.
1294: Logically Collective; No Fortran Support
1296: Input Parameters:
1297: + vec - The vector to modify
1298: . op - The name of the operation
1299: - f - The function that provides the operation.
1301: Level: advanced
1303: Example Usage:
1304: .vb
1305: // some new VecView() implementation, must have the same signature as the function it seeks
1306: // to replace
1307: PetscErrorCode UserVecView(Vec x, PetscViewer viewer)
1308: {
1309: PetscFunctionBeginUser;
1310: // ...
1311: PetscFunctionReturn(PETSC_SUCCESS);
1312: }
1314: // Create a VECMPI which has a pre-defined VecView() implementation
1315: VecCreateMPI(comm, n, N, &x);
1316: // Calls the VECMPI implementation for VecView()
1317: VecView(x, viewer);
1319: VecSetOperation(x, VECOP_VIEW, (PetscErrorCodeFn *)UserVecView);
1320: // Now calls UserVecView()
1321: VecView(x, viewer);
1322: .ve
1324: Notes:
1325: `f` may be `NULL` to remove the operation from `vec`. Depending on the operation this may be
1326: allowed, however some always expect a valid function. In these cases an error will be raised
1327: when calling the interface routine in question.
1329: See `VecOperation` for an up-to-date list of override-able operations. The operations listed
1330: there have the form `VECOP_<OPERATION>`, where `<OPERATION>` is the suffix (in all capital
1331: letters) of the public interface routine (e.g., `VecView()` -> `VECOP_VIEW`).
1333: Overriding a particular `Vec`'s operation has no affect on any other `Vec`s past, present,
1334: or future. The user should also note that overriding a method is "destructive"; the previous
1335: method is not retained in any way.
1337: Each function MUST return `PETSC_SUCCESS` on success and
1338: nonzero on failure.
1340: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecGetOperation()`, `MatSetOperation()`, `MatShellSetOperation()`
1341: @*/
1342: PetscErrorCode VecSetOperation(Vec vec, VecOperation op, PetscErrorCodeFn *f)
1343: {
1344: PetscFunctionBegin;
1346: if (op == VECOP_VIEW && !vec->ops->viewnative) {
1347: vec->ops->viewnative = vec->ops->view;
1348: } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
1349: vec->ops->loadnative = vec->ops->load;
1350: }
1351: ((PetscErrorCodeFn **)vec->ops)[(int)op] = f;
1352: PetscFunctionReturn(PETSC_SUCCESS);
1353: }
1355: /*@
1356: VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1357: used during the assembly process to store values that belong to
1358: other processors.
1360: Not Collective, different processes can have different size stashes
1362: Input Parameters:
1363: + vec - the vector
1364: . size - the initial size of the stash.
1365: - bsize - the initial size of the block-stash(if used).
1367: Options Database Keys:
1368: + -vecstash_initial_size size or size0,size1,...,sizep-1 - set initial size
1369: - -vecstash_block_initial_size bsize or bsize0,bsize1,...,bsizep-1 - set initial block size
1371: Level: intermediate
1373: Notes:
1374: The block-stash is used for values set with `VecSetValuesBlocked()` while
1375: the stash is used for values set with `VecSetValues()`
1377: Run with the option -info and look for output of the form
1378: VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1379: to determine the appropriate value, MM, to use for size and
1380: VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1381: to determine the value, BMM to use for bsize
1383: PETSc attempts to smartly manage the stash size so there is little likelihood setting a
1384: a specific value here will affect performance
1386: .seealso: [](ch_vectors), `Vec`, `VecSetBlockSize()`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecStashView()`
1387: @*/
1388: PetscErrorCode VecStashSetInitialSize(Vec vec, PetscInt size, PetscInt bsize)
1389: {
1390: PetscFunctionBegin;
1392: PetscCall(VecStashSetInitialSize_Private(&vec->stash, size));
1393: PetscCall(VecStashSetInitialSize_Private(&vec->bstash, bsize));
1394: PetscFunctionReturn(PETSC_SUCCESS);
1395: }
1397: /*@
1398: VecSetRandom - Sets all components of a vector to random numbers.
1400: Logically Collective
1402: Input Parameters:
1403: + x - the vector
1404: - rctx - the random number context, formed by `PetscRandomCreate()`, or use `NULL` and it will create one internally.
1406: Output Parameter:
1407: . x - the vector
1409: Example of Usage:
1410: .vb
1411: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1412: VecSetRandom(x,rctx);
1413: PetscRandomDestroy(&rctx);
1414: .ve
1416: Level: intermediate
1418: .seealso: [](ch_vectors), `Vec`, `VecSet()`, `VecSetValues()`, `PetscRandomCreate()`, `PetscRandomDestroy()`
1419: @*/
1420: PetscErrorCode VecSetRandom(Vec x, PetscRandom rctx)
1421: {
1422: PetscRandom randObj = NULL;
1424: PetscFunctionBegin;
1428: VecCheckAssembled(x);
1429: PetscCall(VecSetErrorIfLocked(x, 1));
1431: if (!rctx) {
1432: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)x), &randObj));
1433: PetscCall(PetscRandomSetType(randObj, x->defaultrandtype));
1434: PetscCall(PetscRandomSetFromOptions(randObj));
1435: rctx = randObj;
1436: }
1438: PetscCall(PetscLogEventBegin(VEC_SetRandom, x, rctx, 0, 0));
1439: PetscUseTypeMethod(x, setrandom, rctx);
1440: PetscCall(PetscLogEventEnd(VEC_SetRandom, x, rctx, 0, 0));
1442: PetscCall(PetscRandomDestroy(&randObj));
1443: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1444: PetscFunctionReturn(PETSC_SUCCESS);
1445: }
1447: /*@
1448: VecSetRandomGaussian - Fills a vector with Gaussian random values of the given mean and standard deviation.
1450: Collective
1452: Input Parameters:
1453: + v - the vector to fill
1454: . rng - PETSc random number generator
1455: . mean - desired mean of the Gaussian samples
1456: - std_dev - desired standard deviation
1458: Level: advanced
1460: Note:
1461: For complex builds where `PetscScalar` is complex the imaginary part of all the vector entries is zero
1463: Developer Note:
1464: Uses the Box-Muller transform to generate normally distributed random numbers
1465: from uniform random numbers. Handles edge cases where uniform random values
1466: approach 0 or 1.
1468: .seealso: [](ch_vectors), [](ch_da), `PetscDA`, `PetscRandom`, `PetscRandomSetInterval()`, `VecSetRandom()`
1469: @*/
1470: PetscErrorCode VecSetRandomGaussian(Vec v, PetscRandom rng, PetscReal mean, PetscReal std_dev)
1471: {
1472: PetscInt n, i;
1473: PetscScalar *array;
1474: PetscReal u1, u2;
1475: PetscReal gauss_sample1, gauss_sample2, magnitude, theta;
1476: const PetscReal min_uniform = PETSC_MACHINE_EPSILON;
1477: const PetscInt max_retry_count = 100;
1479: PetscFunctionBegin;
1484: PetscCheck(!PetscIsInfOrNanReal(mean), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Mean must be a finite real number");
1485: PetscCheck(std_dev >= 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Standard deviation must be non-negative, got %g", (double)std_dev);
1486: PetscCheck(!PetscIsInfOrNanReal(std_dev), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Standard deviation must be a finite real number");
1488: PetscCall(VecGetLocalSize(v, &n));
1489: if (n == 0) PetscFunctionReturn(PETSC_SUCCESS);
1491: if (std_dev == 0.0) {
1492: PetscCall(VecSet(v, mean));
1493: PetscFunctionReturn(PETSC_SUCCESS);
1494: }
1496: PetscCall(VecGetArrayWrite(v, &array));
1498: /*
1499: Generate Gaussian-distributed random values using the Box-Muller transform.
1500: This transform converts pairs of uniform random variables U1, U2 ~ Uniform(0,1)
1501: into pairs of independent standard normal variables Z0, Z1 ~ N(0,1):
1502: Z0 = sqrt(-2 * ln(U1)) * cos(2pi * U2)
1503: Z1 = sqrt(-2 * ln(U1)) * sin(2pi * U2)
1504: Then scale and shift to get desired mean and standard deviation.
1505: */
1506: for (i = 0; i < n; i += 2) {
1507: PetscInt retry_count = 0;
1509: /*
1510: Generate U1 and ensure it's not too close to 0 to avoid log(0) singularity.
1511: Add retry limit to prevent infinite loops in case of RNG failure.
1512: */
1513: do {
1514: PetscCall(PetscRandomGetValueReal(rng, &u1));
1515: retry_count++;
1516: PetscCheck(retry_count < max_retry_count, PETSC_COMM_SELF, PETSC_ERR_LIB, "Random number generator failed to produce valid values after %" PetscInt_FMT " attempts", (PetscInt)max_retry_count);
1517: } while (u1 < min_uniform);
1519: PetscCall(PetscRandomGetValueReal(rng, &u2));
1521: /*
1522: Apply Box-Muller transform:
1523: - magnitude: sqrt(-2 * ln(U1)) represents the radial distance from origin
1524: - theta: 2pi * U2 represents the angle uniformly distributed on [0, 2pi]
1525: - Converting from polar to Cartesian coordinates yields two independent samples
1526: */
1527: magnitude = PetscSqrtReal(-2.0 * PetscLogReal(u1));
1528: theta = 2.0 * PETSC_PI * u2;
1529: gauss_sample1 = magnitude * PetscCosReal(theta);
1530: gauss_sample2 = magnitude * PetscSinReal(theta);
1532: /* Scale and shift to achieve desired mean and standard deviation */
1533: array[i] = mean + std_dev * gauss_sample1;
1534: if (i + 1 < n) array[i + 1] = mean + std_dev * gauss_sample2;
1535: }
1537: PetscCall(VecRestoreArrayWrite(v, &array));
1538: PetscFunctionReturn(PETSC_SUCCESS);
1539: }
1541: /*@
1542: VecZeroEntries - puts a `0.0` in each element of a vector
1544: Logically Collective
1546: Input Parameter:
1547: . vec - The vector
1549: Level: beginner
1551: Note:
1552: If the norm of the vector is known to be zero then this skips the unneeded zeroing process
1554: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecSetOptionsPrefix()`, `VecSet()`, `VecSetValues()`
1555: @*/
1556: PetscErrorCode VecZeroEntries(Vec vec)
1557: {
1558: PetscFunctionBegin;
1559: PetscCall(VecSet(vec, 0));
1560: PetscFunctionReturn(PETSC_SUCCESS);
1561: }
1563: /*
1564: VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
1565: processor and a PETSc MPI vector on more than one processor.
1567: Collective
1569: Input Parameter:
1570: . vec - The vector
1572: Level: intermediate
1574: .seealso: [](ch_vectors), `Vec`, `VecSetFromOptions()`, `VecSetType()`
1575: */
1576: static PetscErrorCode VecSetTypeFromOptions_Private(Vec vec, PetscOptionItems PetscOptionsObject)
1577: {
1578: PetscBool opt;
1579: VecType defaultType;
1580: char typeName[256];
1581: PetscMPIInt size;
1583: PetscFunctionBegin;
1584: if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1585: else {
1586: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
1587: if (size > 1) defaultType = VECMPI;
1588: else defaultType = VECSEQ;
1589: }
1591: PetscCall(VecRegisterAll());
1592: PetscCall(PetscOptionsFList("-vec_type", "Vector type", "VecSetType", VecList, defaultType, typeName, 256, &opt));
1593: if (opt) PetscCall(VecSetType(vec, typeName));
1594: else PetscCall(VecSetType(vec, defaultType));
1595: PetscFunctionReturn(PETSC_SUCCESS);
1596: }
1598: /*@
1599: VecSetFromOptions - Configures the vector from the options database.
1601: Collective
1603: Input Parameter:
1604: . vec - The vector
1606: Level: beginner
1608: Notes:
1609: To see all options, run your program with the -help option.
1611: Must be called after `VecCreate()` but before the vector is used.
1613: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecSetOptionsPrefix()`
1614: @*/
1615: PetscErrorCode VecSetFromOptions(Vec vec)
1616: {
1617: PetscBool flg;
1618: PetscInt bind_below = 0;
1620: PetscFunctionBegin;
1623: PetscObjectOptionsBegin((PetscObject)vec);
1624: /* Handle vector type options */
1625: PetscCall(VecSetTypeFromOptions_Private(vec, PetscOptionsObject));
1627: /* Handle specific vector options */
1628: PetscTryTypeMethod(vec, setfromoptions, PetscOptionsObject);
1630: /* Bind to CPU if below a user-specified size threshold.
1631: * This perhaps belongs in the options for the GPU Vec types, but VecBindToCPU() does nothing when called on non-GPU types,
1632: * and putting it here makes is more maintainable than duplicating this for all. */
1633: PetscCall(PetscOptionsInt("-vec_bind_below", "Set the size threshold (in local entries) below which the Vec is bound to the CPU", "VecBindToCPU", bind_below, &bind_below, &flg));
1634: if (flg && vec->map->n < bind_below) PetscCall(VecBindToCPU(vec, PETSC_TRUE));
1636: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1637: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)vec, PetscOptionsObject));
1638: PetscOptionsEnd();
1639: PetscFunctionReturn(PETSC_SUCCESS);
1640: }
1642: /*@
1643: VecSetSizes - Sets the local and global sizes, and checks to determine compatibility of the sizes
1645: Collective
1647: Input Parameters:
1648: + v - the vector
1649: . n - the local size (or `PETSC_DECIDE` to have it set)
1650: - N - the global size (or `PETSC_DETERMINE` to have it set)
1652: Level: intermediate
1654: Notes:
1655: `N` cannot be `PETSC_DETERMINE` if `n` is `PETSC_DECIDE`
1657: If one processor calls this with `N` of `PETSC_DETERMINE` then all processors must, otherwise the program will hang.
1659: If `n` is not `PETSC_DECIDE`, then the value determines the `PetscLayout` of the vector and the ranges returned by
1660: `VecGetOwnershipRange()` and `VecGetOwnershipRanges()`
1662: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecCreateSeq()`, `VecCreateMPI()`, `VecGetSize()`, `PetscSplitOwnership()`, `PetscLayout`,
1663: `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`, `MatSetSizes()`
1664: @*/
1665: PetscErrorCode VecSetSizes(Vec v, PetscInt n, PetscInt N)
1666: {
1667: PetscFunctionBegin;
1669: if (N >= 0) {
1671: PetscCheck(n <= N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local size %" PetscInt_FMT " cannot be larger than global size %" PetscInt_FMT, n, N);
1672: }
1673: PetscCheck(!(v->map->n >= 0 || v->map->N >= 0) || !(v->map->n != n || v->map->N != N), PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot change/reset vector sizes to %" PetscInt_FMT " local %" PetscInt_FMT " global after previously setting them to %" PetscInt_FMT " local %" PetscInt_FMT " global", n, N,
1674: v->map->n, v->map->N);
1675: v->map->n = n;
1676: v->map->N = N;
1677: PetscTryTypeMethod(v, create);
1678: v->ops->create = NULL;
1679: PetscFunctionReturn(PETSC_SUCCESS);
1680: }
1682: /*@
1683: VecSetBlockSize - Sets the block size for future calls to `VecSetValuesBlocked()`
1684: and `VecSetValuesBlockedLocal()`.
1686: Logically Collective
1688: Input Parameters:
1689: + v - the vector
1690: - bs - the blocksize
1692: Level: advanced
1694: Note:
1695: All vectors obtained by `VecDuplicate()` inherit the same blocksize.
1697: Vectors obtained with `DMCreateGlobalVector()` and `DMCreateLocalVector()` generally already have a blocksize set based on the state of the `DM`
1699: .seealso: [](ch_vectors), `Vec`, `VecSetValuesBlocked()`, `VecSetLocalToGlobalMapping()`, `VecGetBlockSize()`
1700: @*/
1701: PetscErrorCode VecSetBlockSize(Vec v, PetscInt bs)
1702: {
1703: PetscFunctionBegin;
1706: PetscCall(PetscLayoutSetBlockSize(v->map, bs));
1707: v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1708: PetscFunctionReturn(PETSC_SUCCESS);
1709: }
1711: /*@
1712: VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for `VecSetValuesBlocked()`
1713: and `VecSetValuesBlockedLocal()`.
1715: Not Collective
1717: Input Parameter:
1718: . v - the vector
1720: Output Parameter:
1721: . bs - the blocksize
1723: Level: advanced
1725: Note:
1726: All vectors obtained by `VecDuplicate()` inherit the same blocksize.
1728: .seealso: [](ch_vectors), `Vec`, `VecSetValuesBlocked()`, `VecSetLocalToGlobalMapping()`, `VecSetBlockSize()`
1729: @*/
1730: PetscErrorCode VecGetBlockSize(Vec v, PetscInt *bs)
1731: {
1732: PetscFunctionBegin;
1734: PetscAssertPointer(bs, 2);
1735: PetscCall(PetscLayoutGetBlockSize(v->map, bs));
1736: PetscFunctionReturn(PETSC_SUCCESS);
1737: }
1739: /*@
1740: VecSetOptionsPrefix - Sets the prefix used for searching for all
1741: `Vec` options in the database.
1743: Logically Collective
1745: Input Parameters:
1746: + v - the `Vec` context
1747: - prefix - the prefix to prepend to all option names
1749: Level: advanced
1751: Note:
1752: A hyphen (-) must NOT be given at the beginning of the prefix name.
1753: The first character of all runtime options is AUTOMATICALLY the hyphen.
1755: .seealso: [](ch_vectors), `Vec`, `VecSetFromOptions()`
1756: @*/
1757: PetscErrorCode VecSetOptionsPrefix(Vec v, const char prefix[])
1758: {
1759: PetscFunctionBegin;
1761: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)v, prefix));
1762: PetscFunctionReturn(PETSC_SUCCESS);
1763: }
1765: /*@
1766: VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1767: `Vec` options in the database.
1769: Logically Collective
1771: Input Parameters:
1772: + v - the `Vec` context
1773: - prefix - the prefix to prepend to all option names
1775: Level: advanced
1777: Note:
1778: A hyphen (-) must NOT be given at the beginning of the prefix name.
1779: The first character of all runtime options is AUTOMATICALLY the hyphen.
1781: .seealso: [](ch_vectors), `Vec`, `VecGetOptionsPrefix()`
1782: @*/
1783: PetscErrorCode VecAppendOptionsPrefix(Vec v, const char prefix[])
1784: {
1785: PetscFunctionBegin;
1787: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)v, prefix));
1788: PetscFunctionReturn(PETSC_SUCCESS);
1789: }
1791: /*@
1792: VecGetOptionsPrefix - Sets the prefix used for searching for all
1793: Vec options in the database.
1795: Not Collective
1797: Input Parameter:
1798: . v - the `Vec` context
1800: Output Parameter:
1801: . prefix - pointer to the prefix string used
1803: Level: advanced
1805: .seealso: [](ch_vectors), `Vec`, `VecAppendOptionsPrefix()`
1806: @*/
1807: PetscErrorCode VecGetOptionsPrefix(Vec v, const char *prefix[])
1808: {
1809: PetscFunctionBegin;
1811: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)v, prefix));
1812: PetscFunctionReturn(PETSC_SUCCESS);
1813: }
1815: /*@C
1816: VecGetState - Gets the state of a `Vec`.
1818: Not Collective
1820: Input Parameter:
1821: . v - the `Vec` context
1823: Output Parameter:
1824: . state - the object state
1826: Level: advanced
1828: Note:
1829: Object state is an integer which gets increased every time
1830: the object is changed. By saving and later querying the object state
1831: one can determine whether information about the object is still current.
1833: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `PetscObjectStateGet()`
1834: @*/
1835: PetscErrorCode VecGetState(Vec v, PetscObjectState *state)
1836: {
1837: PetscFunctionBegin;
1839: PetscAssertPointer(state, 2);
1840: PetscCall(PetscObjectStateGet((PetscObject)v, state));
1841: PetscFunctionReturn(PETSC_SUCCESS);
1842: }
1844: /*@
1845: VecSetUp - Sets up the internal vector data structures for the later use.
1847: Collective
1849: Input Parameter:
1850: . v - the `Vec` context
1852: Level: advanced
1854: Notes:
1855: For basic use of the `Vec` classes the user need not explicitly call
1856: `VecSetUp()`, since these actions will happen automatically.
1858: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecDestroy()`
1859: @*/
1860: PetscErrorCode VecSetUp(Vec v)
1861: {
1862: PetscMPIInt size;
1864: PetscFunctionBegin;
1866: PetscCheck(v->map->n >= 0 || v->map->N >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Sizes not set");
1867: if (!((PetscObject)v)->type_name) {
1868: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1869: if (size == 1) PetscCall(VecSetType(v, VECSEQ));
1870: else PetscCall(VecSetType(v, VECMPI));
1871: }
1872: PetscFunctionReturn(PETSC_SUCCESS);
1873: }
1875: /*
1876: These currently expose the PetscScalar/PetscReal in updating the
1877: cached norm. If we push those down into the implementation these
1878: will become independent of PetscScalar/PetscReal
1879: */
1881: PetscErrorCode VecCopyAsync_Private(Vec x, Vec y, PetscDeviceContext dctx)
1882: {
1883: PetscBool flgs[4];
1884: PetscReal norms[4] = {0.0, 0.0, 0.0, 0.0};
1886: PetscFunctionBegin;
1891: if (x == y) PetscFunctionReturn(PETSC_SUCCESS);
1892: VecCheckSameLocalSize(x, 1, y, 2);
1893: VecCheckAssembled(x);
1894: PetscCall(VecSetErrorIfLocked(y, 2));
1896: #if !defined(PETSC_USE_MIXED_PRECISION)
1897: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
1898: #endif
1900: PetscCall(PetscLogEventBegin(VEC_Copy, x, y, 0, 0));
1901: #if defined(PETSC_USE_MIXED_PRECISION)
1902: extern PetscErrorCode VecGetArray(Vec, double **);
1903: extern PetscErrorCode VecRestoreArray(Vec, double **);
1904: extern PetscErrorCode VecGetArray(Vec, float **);
1905: extern PetscErrorCode VecRestoreArray(Vec, float **);
1906: extern PetscErrorCode VecGetArrayRead(Vec, const double **);
1907: extern PetscErrorCode VecRestoreArrayRead(Vec, const double **);
1908: extern PetscErrorCode VecGetArrayRead(Vec, const float **);
1909: extern PetscErrorCode VecRestoreArrayRead(Vec, const float **);
1910: if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1911: PetscInt i, n;
1912: const float *xx;
1913: double *yy;
1914: PetscCall(VecGetArrayRead(x, &xx));
1915: PetscCall(VecGetArray(y, &yy));
1916: PetscCall(VecGetLocalSize(x, &n));
1917: for (i = 0; i < n; i++) yy[i] = xx[i];
1918: PetscCall(VecRestoreArrayRead(x, &xx));
1919: PetscCall(VecRestoreArray(y, &yy));
1920: } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1921: PetscInt i, n;
1922: float *yy;
1923: const double *xx;
1924: PetscCall(VecGetArrayRead(x, &xx));
1925: PetscCall(VecGetArray(y, &yy));
1926: PetscCall(VecGetLocalSize(x, &n));
1927: for (i = 0; i < n; i++) yy[i] = (float)xx[i];
1928: PetscCall(VecRestoreArrayRead(x, &xx));
1929: PetscCall(VecRestoreArray(y, &yy));
1930: } else PetscUseTypeMethod(x, copy, y);
1931: #else
1932: VecMethodDispatch(x, dctx, VecAsyncFnName(Copy), copy, (Vec, Vec, PetscDeviceContext), y);
1933: #endif
1935: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1936: #if !defined(PETSC_USE_MIXED_PRECISION)
1937: for (PetscInt i = 0; i < 4; i++) {
1938: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)y, NormIds[i], norms[i]));
1939: }
1940: #endif
1942: PetscCall(PetscLogEventEnd(VEC_Copy, x, y, 0, 0));
1943: PetscFunctionReturn(PETSC_SUCCESS);
1944: }
1946: /*@
1947: VecCopy - Copies a vector `y = x`
1949: Logically Collective
1951: Input Parameter:
1952: . x - the vector
1954: Output Parameter:
1955: . y - the copy
1957: Level: beginner
1959: Note:
1960: For default parallel PETSc vectors, both `x` and `y` must be distributed in
1961: the same manner; local copies are done.
1963: Developer Notes:
1964: `PetscCheckSameTypeAndComm`(x,1,y,2) is not used on these vectors because we allow one
1965: of the vectors to be sequential and one to be parallel so long as both have the same
1966: local sizes. This is used in some internal functions in PETSc.
1968: .seealso: [](ch_vectors), `Vec`, `VecDuplicate()`
1969: @*/
1970: PetscErrorCode VecCopy(Vec x, Vec y)
1971: {
1972: PetscFunctionBegin;
1973: PetscCall(VecCopyAsync_Private(x, y, NULL));
1974: PetscFunctionReturn(PETSC_SUCCESS);
1975: }
1977: PetscErrorCode VecSwapAsync_Private(Vec x, Vec y, PetscDeviceContext dctx)
1978: {
1979: PetscReal normxs[4], normys[4];
1980: PetscBool flgxs[4], flgys[4];
1982: PetscFunctionBegin;
1987: PetscCheckSameTypeAndComm(x, 1, y, 2);
1988: VecCheckSameSize(x, 1, y, 2);
1989: VecCheckAssembled(x);
1990: VecCheckAssembled(y);
1991: PetscCall(VecSetErrorIfLocked(x, 1));
1992: PetscCall(VecSetErrorIfLocked(y, 2));
1994: for (PetscInt i = 0; i < 4; i++) {
1995: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], normxs[i], flgxs[i]));
1996: PetscCall(PetscObjectComposedDataGetReal((PetscObject)y, NormIds[i], normys[i], flgys[i]));
1997: }
1999: PetscCall(PetscLogEventBegin(VEC_Swap, x, y, 0, 0));
2000: VecMethodDispatch(x, dctx, VecAsyncFnName(Swap), swap, (Vec, Vec, PetscDeviceContext), y);
2001: PetscCall(PetscLogEventEnd(VEC_Swap, x, y, 0, 0));
2003: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2004: PetscCall(PetscObjectStateIncrease((PetscObject)y));
2005: for (PetscInt i = 0; i < 4; i++) {
2006: if (flgxs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)y, NormIds[i], normxs[i]));
2007: if (flgys[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], normys[i]));
2008: }
2009: PetscFunctionReturn(PETSC_SUCCESS);
2010: }
2011: /*@
2012: VecSwap - Swaps the values between two vectors, `x` and `y`.
2014: Logically Collective
2016: Input Parameters:
2017: + x - the first vector
2018: - y - the second vector
2020: Level: advanced
2022: .seealso: [](ch_vectors), `Vec`, `VecSet()`
2023: @*/
2024: PetscErrorCode VecSwap(Vec x, Vec y)
2025: {
2026: PetscFunctionBegin;
2027: PetscCall(VecSwapAsync_Private(x, y, NULL));
2028: PetscFunctionReturn(PETSC_SUCCESS);
2029: }
2031: /*@
2032: VecStashViewFromOptions - Processes command line options to determine if/how a `VecStash` object is to be viewed.
2034: Collective
2036: Input Parameters:
2037: + obj - the `Vec` containing a stash
2038: . bobj - optional other object that provides the prefix
2039: - name - option to activate viewing
2041: Options Database Key:
2042: . -name [viewertype][:...] - option name and values. See `PetscObjectViewFromOptions()` for the possible arguments
2044: Level: intermediate
2046: Developer Notes:
2047: This cannot use `PetscObjectViewFromOptions()` because it takes a `Vec` as an argument but does not use `VecView()`
2049: .seealso: [](ch_vectors), `Vec`, `VecStashSetInitialSize()`
2050: @*/
2051: PetscErrorCode VecStashViewFromOptions(Vec obj, PetscObject bobj, const char name[])
2052: {
2053: PetscViewer viewer;
2054: PetscBool flg;
2055: PetscViewerFormat format;
2056: char *prefix;
2058: PetscFunctionBegin;
2059: prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
2060: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)obj), ((PetscObject)obj)->options, prefix, name, &viewer, &format, &flg));
2061: if (flg) {
2062: PetscCall(PetscViewerPushFormat(viewer, format));
2063: PetscCall(VecStashView(obj, viewer));
2064: PetscCall(PetscViewerPopFormat(viewer));
2065: PetscCall(PetscViewerDestroy(&viewer));
2066: }
2067: PetscFunctionReturn(PETSC_SUCCESS);
2068: }
2070: /*@
2071: VecStashView - Prints the entries in the vector stash and block stash.
2073: Collective
2075: Input Parameters:
2076: + v - the vector
2077: - viewer - the viewer
2079: Level: advanced
2081: .seealso: [](ch_vectors), `Vec`, `VecSetBlockSize()`, `VecSetValues()`, `VecSetValuesBlocked()`
2082: @*/
2083: PetscErrorCode VecStashView(Vec v, PetscViewer viewer)
2084: {
2085: PetscMPIInt rank;
2086: PetscInt i, j;
2087: PetscBool match;
2088: VecStash *s;
2089: PetscScalar val;
2091: PetscFunctionBegin;
2094: PetscCheckSameComm(v, 1, viewer, 2);
2096: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &match));
2097: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_SUP, "Stash viewer only works with ASCII viewer not %s", ((PetscObject)v)->type_name);
2098: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
2099: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank));
2100: s = &v->bstash;
2102: /* print block stash */
2103: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2104: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Vector Block stash size %" PetscInt_FMT " block size %" PetscInt_FMT "\n", rank, s->n, s->bs));
2105: for (i = 0; i < s->n; i++) {
2106: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " ", rank, s->idx[i]));
2107: for (j = 0; j < s->bs; j++) {
2108: val = s->array[i * s->bs + j];
2109: #if defined(PETSC_USE_COMPLEX)
2110: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "(%18.16e %18.16e) ", (double)PetscRealPart(val), (double)PetscImaginaryPart(val)));
2111: #else
2112: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%18.16e ", (double)val));
2113: #endif
2114: }
2115: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2116: }
2117: PetscCall(PetscViewerFlush(viewer));
2119: s = &v->stash;
2121: /* print basic stash */
2122: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Vector stash size %" PetscInt_FMT "\n", rank, s->n));
2123: for (i = 0; i < s->n; i++) {
2124: val = s->array[i];
2125: #if defined(PETSC_USE_COMPLEX)
2126: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " (%18.16e %18.16e) ", rank, s->idx[i], (double)PetscRealPart(val), (double)PetscImaginaryPart(val)));
2127: #else
2128: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " %18.16e\n", rank, s->idx[i], (double)val));
2129: #endif
2130: }
2131: PetscCall(PetscViewerFlush(viewer));
2132: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2133: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
2134: PetscFunctionReturn(PETSC_SUCCESS);
2135: }
2137: PetscErrorCode PetscOptionsGetVec(PetscOptions options, const char prefix[], const char key[], Vec v, PetscBool *set)
2138: {
2139: PetscInt i, N, rstart, rend;
2140: PetscScalar *xx;
2141: PetscReal *xreal;
2142: PetscBool iset;
2144: PetscFunctionBegin;
2145: PetscCall(VecGetOwnershipRange(v, &rstart, &rend));
2146: PetscCall(VecGetSize(v, &N));
2147: PetscCall(PetscCalloc1(N, &xreal));
2148: PetscCall(PetscOptionsGetRealArray(options, prefix, key, xreal, &N, &iset));
2149: if (iset) {
2150: PetscCall(VecGetArray(v, &xx));
2151: for (i = rstart; i < rend; i++) xx[i - rstart] = xreal[i];
2152: PetscCall(VecRestoreArray(v, &xx));
2153: }
2154: PetscCall(PetscFree(xreal));
2155: if (set) *set = iset;
2156: PetscFunctionReturn(PETSC_SUCCESS);
2157: }
2159: /*@
2160: VecGetLayout - get `PetscLayout` describing a vector layout
2162: Not Collective
2164: Input Parameter:
2165: . x - the vector
2167: Output Parameter:
2168: . map - the layout
2170: Level: developer
2172: Note:
2173: The layout determines what vector elements are contained on each MPI process
2175: .seealso: [](ch_vectors), `PetscLayout`, `Vec`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
2176: @*/
2177: PetscErrorCode VecGetLayout(Vec x, PetscLayout *map)
2178: {
2179: PetscFunctionBegin;
2181: PetscAssertPointer(map, 2);
2182: *map = x->map;
2183: PetscFunctionReturn(PETSC_SUCCESS);
2184: }
2186: /*@
2187: VecSetLayout - set `PetscLayout` describing vector layout
2189: Not Collective
2191: Input Parameters:
2192: + x - the vector
2193: - map - the layout
2195: Level: developer
2197: Note:
2198: It is normally only valid to replace the layout with a layout known to be equivalent.
2200: .seealso: [](ch_vectors), `Vec`, `PetscLayout`, `VecGetLayout()`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
2201: @*/
2202: PetscErrorCode VecSetLayout(Vec x, PetscLayout map)
2203: {
2204: PetscFunctionBegin;
2206: PetscCall(PetscLayoutReference(map, &x->map));
2207: PetscFunctionReturn(PETSC_SUCCESS);
2208: }
2210: /*@
2211: VecFlag - set infinity into the local part of the vector on any subset of MPI processes
2213: Logically Collective
2215: Input Parameters:
2216: + xin - the vector, can be `NULL` but only if on all processes
2217: - flg - indicates if this processes portion of the vector should be set to infinity
2219: Level: developer
2221: Note:
2222: This removes the values from the vector norm cache for all processes by calling `PetscObjectIncrease()`.
2224: This is used for any subset of MPI processes to indicate an failure in a solver, after the next use of `VecNorm()` if
2225: `KSPCheckNorm()` detects an infinity and at least one of the MPI processes has a not converged reason then the `KSP`
2226: object collectively is labeled as not converged.
2228: .seealso: [](ch_vectors), `Vec`, `PetscLayout`, `VecGetLayout()`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
2229: @*/
2230: PetscErrorCode VecFlag(Vec xin, PetscInt flg)
2231: {
2232: // MSVC gives "divide by zero" error at compile time - so declare as volatile to skip this check.
2233: volatile PetscReal one = 1.0, zero = 0.0;
2234: PetscScalar inf;
2236: PetscFunctionBegin;
2237: if (!xin) PetscFunctionReturn(PETSC_SUCCESS);
2239: PetscCall(PetscObjectStateIncrease((PetscObject)xin));
2240: if (flg) {
2241: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2242: inf = one / zero;
2243: PetscCall(PetscFPTrapPop());
2244: if (xin->ops->set) PetscUseTypeMethod(xin, set, inf);
2245: else {
2246: PetscInt n;
2247: PetscScalar *xx;
2249: PetscCall(VecGetLocalSize(xin, &n));
2250: PetscCall(VecGetArrayWrite(xin, &xx));
2251: for (PetscInt i = 0; i < n; ++i) xx[i] = inf;
2252: PetscCall(VecRestoreArrayWrite(xin, &xx));
2253: }
2254: }
2255: PetscFunctionReturn(PETSC_SUCCESS);
2256: }
2258: /*@
2259: VecSetInf - set infinity into the local part of the vector
2261: Not Collective
2263: Input Parameters:
2264: . xin - the vector
2266: Level: developer
2268: Note:
2269: Deprecated, see `VecFlag()`
2270: This is used for any subset of MPI processes to indicate an failure in a solver, after the next use of `VecNorm()` if
2271: `KSPCheckNorm()` detects an infinity and at least one of the MPI processes has a not converged reason then the `KSP`
2272: object collectively is labeled as not converged.
2274: This cannot be called if `xin` has a cached norm available
2276: .seealso: [](ch_vectors), `VecFlag()`, `Vec`, `PetscLayout`, `VecGetLayout()`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
2277: @*/
2278: PetscErrorCode VecSetInf(Vec xin)
2279: {
2280: // MSVC gives "divide by zero" error at compile time - so declare as volatile to skip this check.
2281: volatile PetscReal one = 1.0, zero = 0.0;
2282: PetscScalar inf;
2283: PetscBool flg;
2285: PetscFunctionBegin;
2286: PetscCall(VecNormAvailable(xin, NORM_2, &flg, NULL));
2287: PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot call VecSetInf() if the vector has a cached norm");
2288: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2289: inf = one / zero;
2290: PetscCall(PetscFPTrapPop());
2291: if (xin->ops->set) PetscUseTypeMethod(xin, set, inf);
2292: else {
2293: PetscInt n;
2294: PetscScalar *xx;
2296: PetscCall(VecGetLocalSize(xin, &n));
2297: PetscCall(VecGetArrayWrite(xin, &xx));
2298: for (PetscInt i = 0; i < n; ++i) xx[i] = inf;
2299: PetscCall(VecRestoreArrayWrite(xin, &xx));
2300: }
2301: PetscFunctionReturn(PETSC_SUCCESS);
2302: }
2304: /*@
2305: VecBindToCPU - marks a vector to temporarily stay on the CPU and perform computations on the CPU
2307: Logically collective
2309: Input Parameters:
2310: + v - the vector
2311: - flg - bind to the CPU if value of `PETSC_TRUE`
2313: Level: intermediate
2315: .seealso: [](ch_vectors), `Vec`, `VecBoundToCPU()`
2316: @*/
2317: PetscErrorCode VecBindToCPU(Vec v, PetscBool flg)
2318: {
2319: PetscFunctionBegin;
2322: #if defined(PETSC_HAVE_DEVICE)
2323: if (v->boundtocpu == flg) PetscFunctionReturn(PETSC_SUCCESS);
2324: v->boundtocpu = flg;
2325: PetscTryTypeMethod(v, bindtocpu, flg);
2326: #endif
2327: PetscFunctionReturn(PETSC_SUCCESS);
2328: }
2330: /*@
2331: VecBoundToCPU - query if a vector is bound to the CPU
2333: Not collective
2335: Input Parameter:
2336: . v - the vector
2338: Output Parameter:
2339: . flg - the logical flag
2341: Level: intermediate
2343: .seealso: [](ch_vectors), `Vec`, `VecBindToCPU()`
2344: @*/
2345: PetscErrorCode VecBoundToCPU(Vec v, PetscBool *flg)
2346: {
2347: PetscFunctionBegin;
2349: PetscAssertPointer(flg, 2);
2350: #if defined(PETSC_HAVE_DEVICE)
2351: *flg = v->boundtocpu;
2352: #else
2353: *flg = PETSC_TRUE;
2354: #endif
2355: PetscFunctionReturn(PETSC_SUCCESS);
2356: }
2358: /*@
2359: VecSetBindingPropagates - Sets whether the state of being bound to the CPU for a GPU vector type propagates to child and some other associated objects
2361: Input Parameters:
2362: + v - the vector
2363: - flg - flag indicating whether the boundtocpu flag should be propagated
2365: Level: developer
2367: Notes:
2368: If the value of flg is set to true, then `VecDuplicate()` and `VecDuplicateVecs()` will bind created vectors to GPU if the input vector is bound to the CPU.
2369: The created vectors will also have their bindingpropagates flag set to true.
2371: Developer Notes:
2372: If a `DMDA` has the `-dm_bind_below option` set to true, then vectors created by `DMCreateGlobalVector()` will have `VecSetBindingPropagates()` called on them to
2373: set their bindingpropagates flag to true.
2375: .seealso: [](ch_vectors), `Vec`, `MatSetBindingPropagates()`, `VecGetBindingPropagates()`
2376: @*/
2377: PetscErrorCode VecSetBindingPropagates(Vec v, PetscBool flg)
2378: {
2379: PetscFunctionBegin;
2381: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2382: v->bindingpropagates = flg;
2383: #endif
2384: PetscFunctionReturn(PETSC_SUCCESS);
2385: }
2387: /*@
2388: VecGetBindingPropagates - Gets whether the state of being bound to the CPU for a GPU vector type propagates to child and some other associated objects
2390: Input Parameter:
2391: . v - the vector
2393: Output Parameter:
2394: . flg - flag indicating whether the boundtocpu flag will be propagated
2396: Level: developer
2398: .seealso: [](ch_vectors), `Vec`, `VecSetBindingPropagates()`
2399: @*/
2400: PetscErrorCode VecGetBindingPropagates(Vec v, PetscBool *flg)
2401: {
2402: PetscFunctionBegin;
2404: PetscAssertPointer(flg, 2);
2405: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2406: *flg = v->bindingpropagates;
2407: #else
2408: *flg = PETSC_FALSE;
2409: #endif
2410: PetscFunctionReturn(PETSC_SUCCESS);
2411: }
2413: /*@C
2414: VecSetPinnedMemoryMin - Set the minimum data size for which pinned memory will be used for host (CPU) allocations.
2416: Logically Collective
2418: Input Parameters:
2419: + v - the vector
2420: - mbytes - minimum data size in bytes
2422: Options Database Key:
2423: . -vec_pinned_memory_min size - minimum size (in bytes) for an allocation to use pinned memory on host.
2425: Level: developer
2427: Note:
2428: Specifying -1 ensures that pinned memory will never be used.
2430: .seealso: [](ch_vectors), `Vec`, `VecGetPinnedMemoryMin()`
2431: @*/
2432: PetscErrorCode VecSetPinnedMemoryMin(Vec v, size_t mbytes)
2433: {
2434: PetscFunctionBegin;
2436: #if PetscDefined(HAVE_DEVICE)
2437: v->minimum_bytes_pinned_memory = mbytes;
2438: #endif
2439: PetscFunctionReturn(PETSC_SUCCESS);
2440: }
2442: /*@C
2443: VecGetPinnedMemoryMin - Get the minimum data size for which pinned memory will be used for host (CPU) allocations.
2445: Logically Collective
2447: Input Parameter:
2448: . v - the vector
2450: Output Parameter:
2451: . mbytes - minimum data size in bytes
2453: Level: developer
2455: .seealso: [](ch_vectors), `Vec`, `VecSetPinnedMemoryMin()`
2456: @*/
2457: PetscErrorCode VecGetPinnedMemoryMin(Vec v, size_t *mbytes)
2458: {
2459: PetscFunctionBegin;
2461: PetscAssertPointer(mbytes, 2);
2462: #if PetscDefined(HAVE_DEVICE)
2463: *mbytes = v->minimum_bytes_pinned_memory;
2464: #endif
2465: PetscFunctionReturn(PETSC_SUCCESS);
2466: }
2468: /*@
2469: VecGetOffloadMask - Get the offload mask of a `Vec`
2471: Not Collective
2473: Input Parameter:
2474: . v - the vector
2476: Output Parameter:
2477: . mask - corresponding `PetscOffloadMask` enum value.
2479: Level: intermediate
2481: .seealso: [](ch_vectors), `Vec`, `VecCreateSeqCUDA()`, `VecCreateSeqViennaCL()`, `VecGetArray()`, `VecGetType()`
2482: @*/
2483: PetscErrorCode VecGetOffloadMask(Vec v, PetscOffloadMask *mask)
2484: {
2485: PetscFunctionBegin;
2487: PetscAssertPointer(mask, 2);
2488: *mask = v->offloadmask;
2489: PetscFunctionReturn(PETSC_SUCCESS);
2490: }
2492: #if !defined(PETSC_HAVE_VIENNACL)
2493: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLContext(Vec v, PETSC_UINTPTR_T *ctx)
2494: {
2495: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_context");
2496: }
2498: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLQueue(Vec v, PETSC_UINTPTR_T *queue)
2499: {
2500: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_command_queue");
2501: }
2503: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMem(Vec v, PETSC_UINTPTR_T *queue)
2504: {
2505: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2506: }
2508: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMemRead(Vec v, PETSC_UINTPTR_T *queue)
2509: {
2510: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2511: }
2513: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMemWrite(Vec v, PETSC_UINTPTR_T *queue)
2514: {
2515: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2516: }
2518: PETSC_EXTERN PetscErrorCode VecViennaCLRestoreCLMemWrite(Vec v)
2519: {
2520: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to restore a Vec's cl_mem");
2521: }
2522: #endif
2524: static PetscErrorCode VecErrorWeightedNorms_Basic(Vec U, Vec Y, Vec E, NormType wnormtype, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol, PetscReal ignore_max, PetscReal *norm, PetscInt *norm_loc, PetscReal *norma, PetscInt *norma_loc, PetscReal *normr, PetscInt *normr_loc)
2525: {
2526: const PetscScalar *u, *y;
2527: const PetscScalar *atola = NULL, *rtola = NULL, *erra = NULL;
2528: PetscInt n, n_loc = 0, na_loc = 0, nr_loc = 0;
2529: PetscReal nrm = 0, nrma = 0, nrmr = 0, err_loc[6];
2531: PetscFunctionBegin;
2532: #define SkipSmallValue(a, b, tol) \
2533: if (PetscAbsScalar(a) < tol || PetscAbsScalar(b) < tol) continue
2535: PetscCall(VecGetLocalSize(U, &n));
2536: PetscCall(VecGetArrayRead(U, &u));
2537: PetscCall(VecGetArrayRead(Y, &y));
2538: if (E) PetscCall(VecGetArrayRead(E, &erra));
2539: if (vatol) PetscCall(VecGetArrayRead(vatol, &atola));
2540: if (vrtol) PetscCall(VecGetArrayRead(vrtol, &rtola));
2541: for (PetscInt i = 0; i < n; i++) {
2542: PetscReal err, tol, tola, tolr;
2544: SkipSmallValue(y[i], u[i], ignore_max);
2545: atol = atola ? PetscRealPart(atola[i]) : atol;
2546: rtol = rtola ? PetscRealPart(rtola[i]) : rtol;
2547: err = erra ? PetscAbsScalar(erra[i]) : PetscAbsScalar(y[i] - u[i]);
2548: tola = atol;
2549: tolr = rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
2550: tol = tola + tolr;
2551: if (tola > 0.) {
2552: if (wnormtype == NORM_INFINITY) nrma = PetscMax(nrma, err / tola);
2553: else nrma += PetscSqr(err / tola);
2554: na_loc++;
2555: }
2556: if (tolr > 0.) {
2557: if (wnormtype == NORM_INFINITY) nrmr = PetscMax(nrmr, err / tolr);
2558: else nrmr += PetscSqr(err / tolr);
2559: nr_loc++;
2560: }
2561: if (tol > 0.) {
2562: if (wnormtype == NORM_INFINITY) nrm = PetscMax(nrm, err / tol);
2563: else nrm += PetscSqr(err / tol);
2564: n_loc++;
2565: }
2566: }
2567: if (E) PetscCall(VecRestoreArrayRead(E, &erra));
2568: if (vatol) PetscCall(VecRestoreArrayRead(vatol, &atola));
2569: if (vrtol) PetscCall(VecRestoreArrayRead(vrtol, &rtola));
2570: PetscCall(VecRestoreArrayRead(U, &u));
2571: PetscCall(VecRestoreArrayRead(Y, &y));
2572: #undef SkipSmallValue
2574: err_loc[0] = nrm;
2575: err_loc[1] = nrma;
2576: err_loc[2] = nrmr;
2577: err_loc[3] = (PetscReal)n_loc;
2578: err_loc[4] = (PetscReal)na_loc;
2579: err_loc[5] = (PetscReal)nr_loc;
2580: if (wnormtype == NORM_2) {
2581: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, err_loc, 6, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)U)));
2582: } else {
2583: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, err_loc, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)U)));
2584: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, err_loc + 3, 3, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)U)));
2585: }
2586: if (wnormtype == NORM_2) {
2587: *norm = PetscSqrtReal(err_loc[0]);
2588: *norma = PetscSqrtReal(err_loc[1]);
2589: *normr = PetscSqrtReal(err_loc[2]);
2590: } else {
2591: *norm = err_loc[0];
2592: *norma = err_loc[1];
2593: *normr = err_loc[2];
2594: }
2595: *norm_loc = (PetscInt)err_loc[3];
2596: *norma_loc = (PetscInt)err_loc[4];
2597: *normr_loc = (PetscInt)err_loc[5];
2598: PetscFunctionReturn(PETSC_SUCCESS);
2599: }
2601: /*@
2602: VecErrorWeightedNorms - compute a weighted norm of the difference between two vectors
2604: Collective
2606: Input Parameters:
2607: + U - first vector to be compared
2608: . Y - second vector to be compared
2609: . E - optional third vector representing the error (if not provided, the error is ||U-Y||)
2610: . wnormtype - norm type
2611: . atol - scalar for absolute tolerance
2612: . vatol - vector representing per-entry absolute tolerances (can be ``NULL``)
2613: . rtol - scalar for relative tolerance
2614: . vrtol - vector representing per-entry relative tolerances (can be ``NULL``)
2615: - ignore_max - ignore values smaller than this value in absolute terms.
2617: Output Parameters:
2618: + norm - weighted norm
2619: . norm_loc - number of vector locations used for the weighted norm
2620: . norma - weighted norm based on the absolute tolerance
2621: . norma_loc - number of vector locations used for the absolute weighted norm
2622: . normr - weighted norm based on the relative tolerance
2623: - normr_loc - number of vector locations used for the relative weighted norm
2625: Level: developer
2627: Notes:
2628: This is primarily used for computing weighted local truncation errors in ``TS``.
2630: .seealso: [](ch_vectors), `Vec`, `NormType`, `TSErrorWeightedNorm()`, `TSErrorWeightedENorm()`
2631: @*/
2632: PetscErrorCode VecErrorWeightedNorms(Vec U, Vec Y, Vec E, NormType wnormtype, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol, PetscReal ignore_max, PetscReal *norm, PetscInt *norm_loc, PetscReal *norma, PetscInt *norma_loc, PetscReal *normr, PetscInt *normr_loc)
2633: {
2634: PetscFunctionBegin;
2639: if (E) {
2642: }
2645: if (vatol) {
2648: }
2650: if (vrtol) {
2653: }
2655: PetscAssertPointer(norm, 10);
2656: PetscAssertPointer(norm_loc, 11);
2657: PetscAssertPointer(norma, 12);
2658: PetscAssertPointer(norma_loc, 13);
2659: PetscAssertPointer(normr, 14);
2660: PetscAssertPointer(normr_loc, 15);
2661: PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)U), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
2663: /* There are potentially 5 vectors involved, some of them may happen to be of different type or bound to cpu.
2664: Here we check that they all implement the same operation and call it if so.
2665: Otherwise, we call the _Basic implementation that always works (provided VecGetArrayRead is implemented). */
2666: PetscBool sameop = (PetscBool)(U->ops->errorwnorm && U->ops->errorwnorm == Y->ops->errorwnorm);
2667: if (sameop && E) sameop = (PetscBool)(U->ops->errorwnorm == E->ops->errorwnorm);
2668: if (sameop && vatol) sameop = (PetscBool)(U->ops->errorwnorm == vatol->ops->errorwnorm);
2669: if (sameop && vrtol) sameop = (PetscBool)(U->ops->errorwnorm == vrtol->ops->errorwnorm);
2670: if (sameop) PetscUseTypeMethod(U, errorwnorm, Y, E, wnormtype, atol, vatol, rtol, vrtol, ignore_max, norm, norm_loc, norma, norma_loc, normr, normr_loc);
2671: else PetscCall(VecErrorWeightedNorms_Basic(U, Y, E, wnormtype, atol, vatol, rtol, vrtol, ignore_max, norm, norm_loc, norma, norma_loc, normr, normr_loc));
2672: PetscFunctionReturn(PETSC_SUCCESS);
2673: }