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(PetscOptionsDeprecatedNoObject(PetscObjectComm((PetscObject)vec), ((PetscObject)vec)->prefix, "-vec_view_stash", "-vec_stash_view", "3.26", NULL));
123: PetscCall(VecStashViewFromOptions(vec, NULL, "-vec_stash_view"));
124: PetscCall(PetscLogEventBegin(VEC_AssemblyBegin, vec, 0, 0, 0));
125: PetscTryTypeMethod(vec, assemblybegin);
126: PetscCall(PetscLogEventEnd(VEC_AssemblyBegin, vec, 0, 0, 0));
127: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: /*@
132: VecAssemblyEnd - Completes assembling the vector. This routine should be called after `VecAssemblyBegin()`.
134: Collective
136: Input Parameter:
137: . vec - the vector
139: Options Database Keys:
140: + -vec_view viewer_specification - Call `VecView()` at the conclusion of `VecAssemblyEnd()`. See `PetscOptionsCreateViewer()` for the values of `viewer_specification`.
141: - -vec_stash_view viewer_specification - Call `VecStashView()` during `VecAssemblyBegin()`. See `PetscOptionsCreateViewer()` for the values of `viewer_specification`.
143: Level: beginner
145: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecSetValues()`, `VecView()`, `VecStashView()`, `VecViewFromOptions()`, `VecStashViewFromOptions()`,
146: `PetscObjectViewFromOptions()`
147: @*/
148: PetscErrorCode VecAssemblyEnd(Vec vec)
149: {
150: PetscFunctionBegin;
152: PetscCall(PetscLogEventBegin(VEC_AssemblyEnd, vec, 0, 0, 0));
154: PetscTryTypeMethod(vec, assemblyend);
155: PetscCall(PetscLogEventEnd(VEC_AssemblyEnd, vec, 0, 0, 0));
156: PetscCall(VecViewFromOptions(vec, NULL, "-vec_view"));
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: /*@
161: VecSetPreallocationCOO - set preallocation for a vector using a coordinate format of the entries with global indices
163: Collective
165: Input Parameters:
166: + x - vector being preallocated
167: . ncoo - number of entries
168: - coo_i - entry indices
170: Level: beginner
172: Notes:
173: This and `VecSetValuesCOO()` provide an alternative API to using `VecSetValues()` to provide vector values.
175: This API is particularly efficient for use on GPUs.
177: Entries can be repeated, see `VecSetValuesCOO()`. Negative indices are not allowed unless vector option `VEC_IGNORE_NEGATIVE_INDICES` is set,
178: in which case they, along with the corresponding entries in `VecSetValuesCOO()`, are ignored. If vector option `VEC_NO_OFF_PROC_ENTRIES` is set,
179: remote entries are ignored, otherwise, they will be properly added or inserted to the vector.
181: The array coo_i[] may be freed immediately after calling this function.
183: .seealso: [](ch_vectors), `Vec`, `VecSetValuesCOO()`, `VecSetPreallocationCOOLocal()`
184: @*/
185: PetscErrorCode VecSetPreallocationCOO(Vec x, PetscCount ncoo, const PetscInt coo_i[])
186: {
187: PetscFunctionBegin;
190: if (ncoo) PetscAssertPointer(coo_i, 3);
191: PetscCall(PetscLogEventBegin(VEC_SetPreallocateCOO, x, 0, 0, 0));
192: PetscCall(PetscLayoutSetUp(x->map));
193: if (x->ops->setpreallocationcoo) {
194: PetscUseTypeMethod(x, setpreallocationcoo, ncoo, coo_i);
195: } else {
196: PetscInt ncoo_i;
197: IS is_coo_i;
199: PetscCall(PetscIntCast(ncoo, &ncoo_i));
200: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncoo_i, coo_i, PETSC_COPY_VALUES, &is_coo_i));
201: PetscCall(PetscObjectCompose((PetscObject)x, "__PETSc_coo_i", (PetscObject)is_coo_i));
202: PetscCall(ISDestroy(&is_coo_i));
203: }
204: PetscCall(PetscLogEventEnd(VEC_SetPreallocateCOO, x, 0, 0, 0));
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: /*@
209: VecSetPreallocationCOOLocal - set preallocation for vectors using a coordinate format of the entries with local indices
211: Collective
213: Input Parameters:
214: + x - vector being preallocated
215: . ncoo - number of entries
216: - coo_i - row indices (local numbering; may be modified)
218: Level: beginner
220: Notes:
221: This and `VecSetValuesCOO()` provide an alternative API to using `VecSetValuesLocal()` to provide vector values.
223: This API is particularly efficient for use on GPUs.
225: The local indices are translated using the local to global mapping, thus `VecSetLocalToGlobalMapping()` must have been
226: called prior to this function.
228: The indices coo_i may be modified within this function. They might be translated to corresponding global
229: indices, but the caller should not rely on them having any specific value after this function returns. The arrays
230: can be freed or reused immediately after this function returns.
232: Entries can be repeated. Negative indices and remote indices might be allowed. see `VecSetPreallocationCOO()`.
234: .seealso: [](ch_vectors), `Vec`, `VecSetPreallocationCOO()`, `VecSetValuesCOO()`
235: @*/
236: PetscErrorCode VecSetPreallocationCOOLocal(Vec x, PetscCount ncoo, PetscInt coo_i[])
237: {
238: PetscInt ncoo_i;
239: ISLocalToGlobalMapping ltog;
241: PetscFunctionBegin;
244: if (ncoo) PetscAssertPointer(coo_i, 3);
245: PetscCall(PetscIntCast(ncoo, &ncoo_i));
246: PetscCall(PetscLayoutSetUp(x->map));
247: PetscCall(VecGetLocalToGlobalMapping(x, <og));
248: if (ltog) PetscCall(ISLocalToGlobalMappingApply(ltog, ncoo_i, coo_i, coo_i));
249: PetscCall(VecSetPreallocationCOO(x, ncoo, coo_i));
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
253: /*@
254: VecSetValuesCOO - set values at once in a vector preallocated using `VecSetPreallocationCOO()`
256: Collective
258: Input Parameters:
259: + x - vector being set
260: . coo_v - the value array
261: - imode - the insert mode
263: Level: beginner
265: Note:
266: This and `VecSetPreallocationCOO() or ``VecSetPreallocationCOOLocal()` provide an alternative API to using `VecSetValues()` to provide vector values.
268: This API is particularly efficient for use on GPUs.
270: The values must follow the order of the indices prescribed with `VecSetPreallocationCOO()` or `VecSetPreallocationCOOLocal()`.
271: When repeated entries are specified in the COO indices the `coo_v` values are first properly summed, regardless of the value of `imode`.
272: The imode flag indicates if `coo_v` must be added to the current values of the vector (`ADD_VALUES`) or overwritten (`INSERT_VALUES`).
273: `VecAssemblyBegin()` and `VecAssemblyEnd()` do not need to be called after this routine. It automatically handles the assembly process.
275: .seealso: [](ch_vectors), `Vec`, `VecSetPreallocationCOO()`, `VecSetPreallocationCOOLocal()`, `VecSetValues()`
276: @*/
277: PetscErrorCode VecSetValuesCOO(Vec x, const PetscScalar coo_v[], InsertMode imode)
278: {
279: PetscFunctionBegin;
283: PetscCall(PetscLogEventBegin(VEC_SetValuesCOO, x, 0, 0, 0));
284: if (x->ops->setvaluescoo) {
285: PetscUseTypeMethod(x, setvaluescoo, coo_v, imode);
286: PetscCall(PetscObjectStateIncrease((PetscObject)x));
287: } else {
288: IS is_coo_i;
289: const PetscInt *coo_i;
290: PetscInt ncoo;
291: PetscMemType mtype;
293: PetscCall(PetscGetMemType(coo_v, &mtype));
294: PetscCheck(mtype == PETSC_MEMTYPE_HOST, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_WRONG, "The basic VecSetValuesCOO() only supports v[] on host");
295: PetscCall(PetscObjectQuery((PetscObject)x, "__PETSc_coo_i", (PetscObject *)&is_coo_i));
296: PetscCheck(is_coo_i, PetscObjectComm((PetscObject)x), PETSC_ERR_COR, "Missing coo_i IS");
297: PetscCall(ISGetLocalSize(is_coo_i, &ncoo));
298: PetscCall(ISGetIndices(is_coo_i, &coo_i));
299: if (imode != ADD_VALUES) PetscCall(VecZeroEntries(x));
300: PetscCall(VecSetValues(x, ncoo, coo_i, coo_v, ADD_VALUES));
301: PetscCall(ISRestoreIndices(is_coo_i, &coo_i));
302: PetscCall(VecAssemblyBegin(x));
303: PetscCall(VecAssemblyEnd(x));
304: }
305: PetscCall(PetscLogEventEnd(VEC_SetValuesCOO, x, 0, 0, 0));
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: 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))
310: {
311: PetscErrorCode (*async_fn)(Vec, Vec, Vec, PetscDeviceContext) = NULL;
313: PetscFunctionBegin;
320: PetscCheckSameTypeAndComm(x, 2, y, 3);
321: PetscCheckSameTypeAndComm(y, 3, w, 1);
322: VecCheckSameSize(w, 1, x, 2);
323: VecCheckSameSize(w, 1, y, 3);
324: VecCheckAssembled(x);
325: VecCheckAssembled(y);
326: PetscCall(VecSetErrorIfLocked(w, 1));
329: if (dctx) PetscCall(PetscObjectQueryFunction((PetscObject)w, async_name, &async_fn));
330: if (event) PetscCall(PetscLogEventBegin(event, x, y, w, 0));
331: if (async_fn) PetscCall((*async_fn)(w, x, y, dctx));
332: else PetscCall((*pointwise_op)(w, x, y));
333: if (event) PetscCall(PetscLogEventEnd(event, x, y, w, 0));
334: PetscCall(PetscObjectStateIncrease((PetscObject)w));
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: PetscErrorCode VecPointwiseMaxAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
339: {
340: PetscFunctionBegin;
341: // REVIEW ME: no log event?
342: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMax), w->ops->pointwisemax));
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: /*@
347: VecPointwiseMax - Computes the component-wise maximum `w[i] = max(x[i], y[i])`.
349: Logically Collective
351: Input Parameters:
352: + x - the first input vector
353: - y - the second input vector
355: Output Parameter:
356: . w - the result
358: Level: advanced
360: Notes:
361: Any subset of the `x`, `y`, and `w` may be the same vector.
363: For complex numbers compares only the real part
365: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
366: @*/
367: PetscErrorCode VecPointwiseMax(Vec w, Vec x, Vec y)
368: {
369: PetscFunctionBegin;
370: PetscCall(VecPointwiseMaxAsync_Private(w, x, y, NULL));
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: PetscErrorCode VecPointwiseMinAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
375: {
376: PetscFunctionBegin;
377: // REVIEW ME: no log event?
378: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMin), w->ops->pointwisemin));
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: /*@
383: VecPointwiseMin - Computes the component-wise minimum `w[i] = min(x[i], y[i])`.
385: Logically Collective
387: Input Parameters:
388: + x - the first input vector
389: - y - the second input vector
391: Output Parameter:
392: . w - the result
394: Level: advanced
396: Notes:
397: Any subset of the `x`, `y`, and `w` may be the same vector.
399: For complex numbers compares only the real part
401: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
402: @*/
403: PetscErrorCode VecPointwiseMin(Vec w, Vec x, Vec y)
404: {
405: PetscFunctionBegin;
406: PetscCall(VecPointwiseMinAsync_Private(w, x, y, NULL));
407: PetscFunctionReturn(PETSC_SUCCESS);
408: }
410: PetscErrorCode VecPointwiseMaxAbsAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
411: {
412: PetscFunctionBegin;
413: // REVIEW ME: no log event?
414: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMaxAbs), w->ops->pointwisemaxabs));
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: /*@
419: VecPointwiseMaxAbs - Computes the component-wise maximum of the absolute values `w[i] = max(abs(x[i]), abs(y[i]))`.
421: Logically Collective
423: Input Parameters:
424: + x - the first input vector
425: - y - the second input vector
427: Output Parameter:
428: . w - the result
430: Level: advanced
432: Notes:
433: Any subset of the `x`, `y`, and `w` may be the same vector.
435: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMin()`, `VecPointwiseMax()`, `VecMaxPointwiseDivide()`
436: @*/
437: PetscErrorCode VecPointwiseMaxAbs(Vec w, Vec x, Vec y)
438: {
439: PetscFunctionBegin;
440: PetscCall(VecPointwiseMaxAbsAsync_Private(w, x, y, NULL));
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
444: PetscErrorCode VecPointwiseDivideAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
445: {
446: PetscFunctionBegin;
447: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, VEC_PointwiseDivide, VecAsyncFnName(PointwiseDivide), w->ops->pointwisedivide));
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: /*@
452: VecPointwiseDivide - Computes the component-wise division `w[i] = x[i] / y[i]`.
454: Logically Collective
456: Input Parameters:
457: + x - the numerator vector
458: - y - the denominator vector
460: Output Parameter:
461: . w - the result
463: Level: advanced
465: Notes:
466: Any subset of the `x`, `y`, and `w` may be the same vector.
468: If a particular `y[i]` is zero and `x[i]` is also zero, `w[i]` is set to one. If instead `x[i]` is not zero, then `w[i]` is zero.
470: .seealso: [](ch_vectors), `Vec`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
471: @*/
472: PetscErrorCode VecPointwiseDivide(Vec w, Vec x, Vec y)
473: {
474: PetscFunctionBegin;
475: PetscCall(VecPointwiseDivideAsync_Private(w, x, y, NULL));
476: PetscFunctionReturn(PETSC_SUCCESS);
477: }
479: #define VEC_POINTWISE_SIGN_LOOP(y, x, n, func) \
480: PetscPragmaSIMD \
481: for (PetscInt i = 0; i < (n); i++) (y)[i] = func(PetscRealPart((x)[i]))
483: #define VEC_POINTWISE_SIGN_DISPATCH(y, x, n, sign_type) \
484: do { \
485: switch (sign_type) { \
486: case VEC_SIGN_ZERO_TO_ZERO: \
487: VEC_POINTWISE_SIGN_LOOP(y, x, n, VecSignZeroToZero_Private); \
488: break; \
489: case VEC_SIGN_ZERO_TO_SIGNED_ZERO: \
490: VEC_POINTWISE_SIGN_LOOP(y, x, n, VecSignZeroToSignedZero_Private); \
491: break; \
492: case VEC_SIGN_ZERO_TO_SIGNED_UNIT: \
493: VEC_POINTWISE_SIGN_LOOP(y, x, n, VecSignZeroToSignedUnit_Private); \
494: break; \
495: default: \
496: PetscUnreachable(); \
497: } \
498: } while (0)
500: PetscErrorCode VecPointwiseSignAsync_Private(Vec y, Vec x, VecSignMode sign_type, PetscDeviceContext dctx)
501: {
502: PetscOffloadMask mask;
503: PetscBool is_host;
504: PetscErrorCode (*async_fn)(Vec, Vec, VecSignMode, PetscDeviceContext) = NULL;
506: PetscFunctionBegin;
511: VecCheckSameSize(y, 1, x, 2);
512: VecCheckAssembled(x);
513: VecCheckAssembled(y);
514: PetscCall(VecSetErrorIfLocked(y, 1));
516: PetscCall(VecGetOffloadMask(x, &mask));
517: is_host = PetscOffloadHost(mask) ? PETSC_TRUE : PETSC_FALSE;
518: if (!is_host) PetscCall(PetscObjectQueryFunction((PetscObject)y, VEC_ASYNC_FN_NAME("PointwiseSign"), &async_fn));
519: if (async_fn) PetscCall((*async_fn)(y, x, sign_type, dctx));
520: else {
521: PetscInt n;
523: PetscCall(VecGetLocalSize(y, &n));
524: if (y == x) {
525: PetscScalar *_y;
527: PetscCall(VecGetArray(y, &_y));
528: VEC_POINTWISE_SIGN_DISPATCH(_y, _y, n, sign_type);
529: PetscCall(VecRestoreArray(y, &_y));
530: } else {
531: PetscScalar *_y;
532: const PetscScalar *_x;
534: PetscCall(VecGetArrayWrite(y, &_y));
535: PetscCall(VecGetArrayRead(x, &_x));
536: VEC_POINTWISE_SIGN_DISPATCH(_y, _x, n, sign_type);
537: PetscCall(VecRestoreArrayRead(x, &_x));
538: PetscCall(VecRestoreArrayWrite(y, &_y));
539: }
540: }
541: PetscCall(PetscObjectStateIncrease((PetscObject)y));
542: PetscFunctionReturn(PETSC_SUCCESS);
543: }
545: /*@
546: VecPointwiseSign - Computes the component-wise sign `y[i] = sign(x[i])`.
548: Logically Collective
550: Input Parameters:
551: + x - the input vector
552: - sign_type - `VecSignMode` indicating how the function should map zero values.
554: Output Parameter:
555: . y - the sign vector of `x`
557: Level: beginner
559: .seealso: [](ch_vectors), `Vec`, `VecSignMode`
560: @*/
561: PetscErrorCode VecPointwiseSign(Vec y, Vec x, VecSignMode sign_type)
562: {
563: PetscFunctionBegin;
564: PetscCall(VecPointwiseSignAsync_Private(y, x, sign_type, NULL));
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: PetscErrorCode VecPointwiseMultAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
569: {
570: PetscFunctionBegin;
572: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, VEC_PointwiseMult, VecAsyncFnName(PointwiseMult), w->ops->pointwisemult));
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: /*@
577: VecPointwiseMult - Computes the component-wise multiplication `w[i] = x[i] * y[i]`.
579: Logically Collective
581: Input Parameters:
582: + x - the first vector
583: - y - the second vector
585: Output Parameter:
586: . w - the result
588: Level: advanced
590: Note:
591: Any subset of the `x`, `y`, and `w` may be the same vector.
593: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
594: @*/
595: PetscErrorCode VecPointwiseMult(Vec w, Vec x, Vec y)
596: {
597: PetscFunctionBegin;
598: PetscCall(VecPointwiseMultAsync_Private(w, x, y, NULL));
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: /*@
603: VecDuplicate - Creates a new vector of the same type as an existing vector.
605: Collective
607: Input Parameter:
608: . v - a vector to mimic
610: Output Parameter:
611: . newv - location to put new vector
613: Level: beginner
615: Notes:
616: `VecDuplicate()` DOES NOT COPY the vector entries, but rather allocates storage
617: for the new vector. Use `VecCopy()` to copy a vector.
619: PETSc `Vec` always have all zero entries when created with `VecDuplicate()` until routines such as `VecSet()` or `VecSetValues()`
620: are used to change the values. There is no reason to call `VecZeroEntries()` after creation.
622: Use `VecDestroy()` to free the space. Use `VecDuplicateVecs()` to get several
623: vectors.
625: .seealso: [](ch_vectors), `Vec`, `VecDestroy()`, `VecDuplicateVecs()`, `VecCreate()`, `VecCopy()`
626: @*/
627: PetscErrorCode VecDuplicate(Vec v, Vec *newv)
628: {
629: PetscFunctionBegin;
631: PetscAssertPointer(newv, 2);
633: PetscUseTypeMethod(v, duplicate, newv);
634: #if PetscDefined(HAVE_DEVICE)
635: if (v->boundtocpu && v->bindingpropagates) {
636: PetscCall(VecSetBindingPropagates(*newv, PETSC_TRUE));
637: PetscCall(VecBindToCPU(*newv, PETSC_TRUE));
638: }
639: #endif
640: PetscCall(PetscObjectStateIncrease((PetscObject)*newv));
641: PetscFunctionReturn(PETSC_SUCCESS);
642: }
644: /*@
645: VecDestroy - Destroys a vector.
647: Collective
649: Input Parameter:
650: . v - the vector
652: Level: beginner
654: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecDuplicate()`, `VecDestroyVecs()`
655: @*/
656: PetscErrorCode VecDestroy(Vec *v)
657: {
658: PetscFunctionBegin;
659: PetscAssertPointer(v, 1);
660: if (!*v) PetscFunctionReturn(PETSC_SUCCESS);
662: if (--((PetscObject)*v)->refct > 0) {
663: *v = NULL;
664: PetscFunctionReturn(PETSC_SUCCESS);
665: }
667: PetscCall(PetscObjectSAWsViewOff((PetscObject)*v));
668: /* destroy the internal part */
669: PetscTryTypeMethod(*v, destroy);
670: PetscCall(PetscFree((*v)->defaultrandtype));
671: /* destroy the external/common part */
672: PetscCall(PetscLayoutDestroy(&(*v)->map));
673: PetscCall(PetscHeaderDestroy(v));
674: PetscFunctionReturn(PETSC_SUCCESS);
675: }
677: /*@
678: VecDuplicateVecs - Creates several vectors of the same type as an existing vector.
680: Collective
682: Input Parameters:
683: + m - the number of vectors to obtain
684: - v - a vector to mimic
686: Output Parameter:
687: . V - location to put pointer to array of vectors
689: Level: intermediate
691: Notes:
692: Use `VecDestroyVecs()` to free the space. Use `VecDuplicate()` to form a single
693: vector.
695: PETSc `Vec` always have all zero entries when created with `VecDuplicateVecs()` until routines such as `VecSet()` or `VecSetValues()`
696: are used to change the values. There is no reason to call `VecZeroEntries()` after creation.
698: Some implementations ensure that the arrays accessed by each vector are contiguous in memory. Certain `VecMDot()` and `VecMAXPY()`
699: implementations utilize this property to use BLAS 2 operations for higher efficiency. This is especially useful in `KSPGMRES`, see
700: `KSPGMRESSetPreAllocateVectors()`.
702: Fortran Note:
703: .vb
704: Vec, pointer :: V(:)
705: .ve
707: .seealso: [](ch_vectors), `Vec`, [](ch_fortran), `VecDestroyVecs()`, `VecDuplicate()`, `VecCreate()`, `VecMDot()`, `VecMAXPY()`, `KSPGMRES`,
708: `KSPGMRESSetPreAllocateVectors()`
709: @*/
710: PetscErrorCode VecDuplicateVecs(Vec v, PetscInt m, Vec *V[])
711: {
712: PetscFunctionBegin;
714: PetscAssertPointer(V, 3);
716: PetscUseTypeMethod(v, duplicatevecs, m, V);
717: #if PetscDefined(HAVE_VIENNACL) || PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
718: if (v->boundtocpu && v->bindingpropagates) {
719: for (PetscInt i = 0; i < m; i++) {
720: /* Since ops->duplicatevecs might itself propagate the value of boundtocpu,
721: * avoid unnecessary overhead by only calling VecBindToCPU() if the vector isn't already bound. */
722: if (!(*V)[i]->boundtocpu) {
723: PetscCall(VecSetBindingPropagates((*V)[i], PETSC_TRUE));
724: PetscCall(VecBindToCPU((*V)[i], PETSC_TRUE));
725: }
726: }
727: }
728: #endif
729: PetscFunctionReturn(PETSC_SUCCESS);
730: }
732: /*@
733: VecDestroyVecs - Frees a block of vectors obtained with `VecDuplicateVecs()`.
735: Collective
737: Input Parameters:
738: + m - the number of vectors previously obtained, if zero no vectors are destroyed
739: - vv - pointer to pointer to array of vector pointers, if `NULL` no vectors are destroyed
741: Level: intermediate
743: .seealso: [](ch_vectors), `Vec`, [](ch_fortran), `VecDuplicateVecs()`, `VecDestroyVecsf90()`
744: @*/
745: PetscErrorCode VecDestroyVecs(PetscInt m, Vec *vv[])
746: {
747: PetscFunctionBegin;
748: PetscAssertPointer(vv, 2);
749: PetscCheck(m >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Trying to destroy negative number of vectors %" PetscInt_FMT, m);
750: if (!m || !*vv) {
751: *vv = NULL;
752: PetscFunctionReturn(PETSC_SUCCESS);
753: }
756: PetscCall((*(**vv)->ops->destroyvecs)(m, *vv));
757: *vv = NULL;
758: PetscFunctionReturn(PETSC_SUCCESS);
759: }
761: /*@
762: VecViewFromOptions - View a vector based on values in the options database
764: Collective
766: Input Parameters:
767: + A - the vector
768: . obj - optional object that provides the options prefix for this viewing, use `NULL` to use the prefix of `A`
769: - name - command line option
771: Options Database Key:
772: . -name viewer_specification - See `PetscOptionsCreateViewer()` for the values of `viewer_specification`
774: Level: intermediate
776: Note:
777: This checks the options database, creates the viewer on-the-fly, uses it and then destroys it. Hence it should not be called in heavily used routines,
778: rather `PetscOptionsCreateViewer()` should be used to construct the viewer once which can then be utilized in the heavily used routine.
780: .seealso: [](ch_vectors), `Vec`, `VecView()`, `PetscObjectViewFromOptions()`, `PetscOptionsCreateViewer()`, `VecCreate()`
781: @*/
782: PetscErrorCode VecViewFromOptions(Vec A, PeOp PetscObject obj, const char name[])
783: {
784: PetscFunctionBegin;
786: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
787: PetscFunctionReturn(PETSC_SUCCESS);
788: }
790: /*@
791: VecView - Views a vector object.
793: Collective
795: Input Parameters:
796: + vec - the vector
797: - viewer - an optional `PetscViewer` visualization context
799: Options Database Key:
800: . -vec_view viewer_specification - Call `VecView()` at the conclusion of `VecAssemblyEnd()`. See `PetscOptionsCreateViewer()` for the values of `viewer_specification`.
802: Level: beginner
804: Notes:
805: The available visualization contexts include
806: + `PETSC_VIEWER_STDOUT_SELF` - for sequential vectors
807: . `PETSC_VIEWER_STDOUT_WORLD` - for parallel vectors created on `PETSC_COMM_WORLD`
808: - `PETSC_VIEWER_STDOUT`_(comm) - for parallel vectors created on MPI communicator comm
810: You can change the format the vector is printed using the
811: option `PetscViewerPushFormat()`.
813: The user can open alternative viewers with
814: + `PetscViewerASCIIOpen()` - Outputs vector to a specified file
815: . `PetscViewerBinaryOpen()` - Outputs vector in binary to a
816: specified file; corresponding input uses `VecLoad()`
817: . `PetscViewerDrawOpen()` - Outputs vector to an X window display
818: . `PetscViewerSocketOpen()` - Outputs vector to Socket viewer
819: - `PetscViewerHDF5Open()` - Outputs vector to HDF5 file viewer
821: The user can call `PetscViewerPushFormat()` to specify the output
822: format of ASCII printed objects (when using `PETSC_VIEWER_STDOUT_SELF`,
823: `PETSC_VIEWER_STDOUT_WORLD` and `PetscViewerASCIIOpen()`). Available formats include
824: + `PETSC_VIEWER_DEFAULT` - default, prints vector contents
825: . `PETSC_VIEWER_ASCII_MATLAB` - prints vector contents in MATLAB format
826: . `PETSC_VIEWER_ASCII_INDEX` - prints vector contents, including indices of vector elements
827: - `PETSC_VIEWER_ASCII_COMMON` - prints vector contents, using a
828: format common among all vector types
830: `VecViewFromOptions()` provides an alternative to this routine that only views the vector if the requested value
831: is provided in the options database.
833: You can pass any number of vector objects, or other PETSc objects to the same viewer.
835: In the debugger you can do call `VecView`(v,0) to display the vector. (The same holds for any PETSc object viewer).
837: Notes for binary viewer:
838: If you pass multiple vectors to a binary viewer you can read them back in the same order
839: with `VecLoad()`.
841: If the blocksize of the vector is greater than one then you must provide a unique prefix to
842: the vector with `PetscObjectSetOptionsPrefix`((`PetscObject`)vec,"uniqueprefix"); BEFORE calling `VecView()` on the
843: vector to be stored and then set that same unique prefix on the vector that you pass to `VecLoad()`. The blocksize
844: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
845: filename. If you copy the binary file, make sure you copy the associated .info file with it.
847: See the manual page for `VecLoad()` on the exact format the binary viewer stores
848: the values in the file.
850: Notes for HDF5 Viewer:
851: The name of the `Vec` (given with `PetscObjectSetName()` is the name that is used
852: for the object in the HDF5 file. If you wish to store the same Vec into multiple
853: datasets in the same file (typically with different values), you must change its
854: name each time before calling the `VecView()`. To load the same vector,
855: the name of the Vec object passed to `VecLoad()` must be the same.
857: If the block size of the vector is greater than 1 then it is used as the first dimension in the HDF5 array.
858: If the function `PetscViewerHDF5SetBaseDimension2()`is called then even if the block size is one it will
859: be used as the first dimension in the HDF5 array (that is the HDF5 array will always be two dimensional)
860: See also `PetscViewerHDF5SetTimestep()` which adds an additional complication to reading and writing `Vec`
861: with the HDF5 viewer.
863: .seealso: [](ch_vectors), `Vec`, `VecViewFromOptions()`, `PetscViewerASCIIOpen()`, `PetscViewerDrawOpen()`, `PetscDrawLGCreate()`,
864: `PetscViewerSocketOpen()`, `PetscViewerBinaryOpen()`, `VecLoad()`, `PetscViewerCreate()`,
865: `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`, `PetscViewerHDF5SetTimestep()`, `PetscOptionsCreateViewer()`
866: @*/
867: PetscErrorCode VecView(Vec vec, PetscViewer viewer)
868: {
869: PetscBool isascii;
870: PetscViewerFormat format;
871: PetscMPIInt size;
873: PetscFunctionBegin;
876: VecCheckAssembled(vec);
877: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec), &viewer));
879: PetscCall(PetscViewerGetFormat(viewer, &format));
880: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
881: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) PetscFunctionReturn(PETSC_SUCCESS);
883: PetscCheck(!vec->stash.n && !vec->bstash.n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call VecAssemblyBegin/End() before viewing this vector");
885: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
886: if (isascii) {
887: PetscInt rows, bs;
889: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)vec, viewer));
890: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
891: PetscCall(PetscViewerASCIIPushTab(viewer));
892: PetscCall(VecGetSize(vec, &rows));
893: PetscCall(VecGetBlockSize(vec, &bs));
894: if (bs != 1) {
895: PetscCall(PetscViewerASCIIPrintf(viewer, "length=%" PetscInt_FMT ", bs=%" PetscInt_FMT "\n", rows, bs));
896: } else {
897: PetscCall(PetscViewerASCIIPrintf(viewer, "length=%" PetscInt_FMT "\n", rows));
898: }
899: PetscCall(PetscViewerASCIIPopTab(viewer));
900: }
901: }
902: PetscCall(VecLockReadPush(vec));
903: PetscCall(PetscLogEventBegin(VEC_View, vec, viewer, 0, 0));
904: if ((format == PETSC_VIEWER_NATIVE || format == PETSC_VIEWER_LOAD_BALANCE) && vec->ops->viewnative) {
905: PetscUseTypeMethod(vec, viewnative, viewer);
906: } else {
907: PetscUseTypeMethod(vec, view, viewer);
908: }
909: PetscCall(VecLockReadPop(vec));
910: PetscCall(PetscLogEventEnd(VEC_View, vec, viewer, 0, 0));
911: PetscFunctionReturn(PETSC_SUCCESS);
912: }
914: #if PetscDefined(USE_DEBUG)
915: #include <../src/sys/totalview/tv_data_display.h>
916: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
917: {
918: const PetscScalar *values;
919: char type[32];
921: TV_add_row("Local rows", "int", &v->map->n);
922: TV_add_row("Global rows", "int", &v->map->N);
923: TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
924: PetscCall(VecGetArrayRead((Vec)v, &values));
925: PetscCall(PetscSNPrintf(type, 32, "double[%" PetscInt_FMT "]", v->map->n));
926: TV_add_row("values", type, values);
927: PetscCall(VecRestoreArrayRead((Vec)v, &values));
928: return TV_format_OK;
929: }
930: #endif
932: /*@
933: VecViewNative - Views a vector object with the original type specific viewer
935: Collective
937: Input Parameters:
938: + vec - the vector
939: - viewer - an optional `PetscViewer` visualization context
941: Level: developer
943: Note:
944: This can be used with, for example, vectors obtained with `DMCreateGlobalVector()` for a `DMDA` to display the vector
945: in the PETSc storage format (each MPI process values follow the previous MPI processes) instead of the "natural" grid
946: ordering.
948: .seealso: [](ch_vectors), `Vec`, `PetscViewerASCIIOpen()`, `PetscViewerDrawOpen()`, `PetscDrawLGCreate()`, `VecView()`,
949: `PetscViewerSocketOpen()`, `PetscViewerBinaryOpen()`, `VecLoad()`, `PetscViewerCreate()`,
950: `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`, `PetscViewerHDF5SetTimestep()`
951: @*/
952: PetscErrorCode VecViewNative(Vec vec, PetscViewer viewer)
953: {
954: PetscFunctionBegin;
957: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec), &viewer));
959: PetscUseTypeMethod(vec, viewnative, viewer);
960: PetscFunctionReturn(PETSC_SUCCESS);
961: }
963: /*@
964: VecGetSize - Returns the global number of elements of the vector.
966: Not Collective
968: Input Parameter:
969: . x - the vector
971: Output Parameter:
972: . size - the global length of the vector
974: Level: beginner
976: .seealso: [](ch_vectors), `Vec`, `VecGetLocalSize()`
977: @*/
978: PetscErrorCode VecGetSize(Vec x, PetscInt *size)
979: {
980: PetscFunctionBegin;
982: PetscAssertPointer(size, 2);
984: PetscUseTypeMethod(x, getsize, size);
985: PetscFunctionReturn(PETSC_SUCCESS);
986: }
988: /*@
989: VecGetLocalSize - Returns the number of elements of the vector stored
990: in local memory (that is on this MPI process)
992: Not Collective
994: Input Parameter:
995: . x - the vector
997: Output Parameter:
998: . size - the length of the local piece of the vector
1000: Level: beginner
1002: .seealso: [](ch_vectors), `Vec`, `VecGetSize()`
1003: @*/
1004: PetscErrorCode VecGetLocalSize(Vec x, PetscInt *size)
1005: {
1006: PetscFunctionBegin;
1008: PetscAssertPointer(size, 2);
1010: PetscUseTypeMethod(x, getlocalsize, size);
1011: PetscFunctionReturn(PETSC_SUCCESS);
1012: }
1014: /*@
1015: VecGetOwnershipRange - Returns the range of indices owned by
1016: this process. The vector is laid out with the
1017: first `n1` elements on the first processor, next `n2` elements on the
1018: second, etc. For certain parallel layouts this range may not be
1019: well defined.
1021: Not Collective
1023: Input Parameter:
1024: . x - the vector
1026: Output Parameters:
1027: + low - the first local element, pass in `NULL` if not interested
1028: - high - one more than the last local element, pass in `NULL` if not interested
1030: Level: beginner
1032: Notes:
1033: If the `Vec` was obtained from a `DM` with `DMCreateGlobalVector()`, then the range values are determined by the specific `DM`.
1035: If the `Vec` was created directly the range values are determined by the local size passed to `VecSetSizes()` or `VecCreateMPI()`.
1036: If `PETSC_DECIDE` was passed as the local size, then the vector uses default values for the range using `PetscSplitOwnership()`.
1038: The high argument is one more than the last element stored locally.
1040: For certain `DM`, such as `DMDA`, it is better to use `DM` specific routines, such as `DMDAGetGhostCorners()`, to determine
1041: the local values in the vector.
1043: .seealso: [](ch_vectors), `Vec`, `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `VecGetOwnershipRanges()`, `PetscSplitOwnership()`,
1044: `VecSetSizes()`, `VecCreateMPI()`, `PetscLayout`, `DMDAGetGhostCorners()`, `DM`
1045: @*/
1046: PetscErrorCode VecGetOwnershipRange(Vec x, PetscInt *low, PetscInt *high)
1047: {
1048: PetscFunctionBegin;
1051: if (low) PetscAssertPointer(low, 2);
1052: if (high) PetscAssertPointer(high, 3);
1053: if (low) *low = x->map->rstart;
1054: if (high) *high = x->map->rend;
1055: PetscFunctionReturn(PETSC_SUCCESS);
1056: }
1058: /*@
1059: VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
1060: The vector is laid out with the
1061: first `n1` elements on the first processor, next `n2` elements on the
1062: second, etc. For certain parallel layouts this range may not be
1063: well defined.
1065: Not Collective
1067: Input Parameter:
1068: . x - the vector
1070: Output Parameter:
1071: . ranges - array of length `size` + 1 with the start and end+1 for each process
1073: Level: beginner
1075: Notes:
1076: If the `Vec` was obtained from a `DM` with `DMCreateGlobalVector()`, then the range values are determined by the specific `DM`.
1078: If the `Vec` was created directly the range values are determined by the local size passed to `VecSetSizes()` or `VecCreateMPI()`.
1079: If `PETSC_DECIDE` was passed as the local size, then the vector uses default values for the range using `PetscSplitOwnership()`.
1081: The high argument is one more than the last element stored locally.
1083: For certain `DM`, such as `DMDA`, it is better to use `DM` specific routines, such as `DMDAGetGhostCorners()`, to determine
1084: the local values in the vector.
1086: The high argument is one more than the last element stored locally.
1088: If `ranges` are used after all vectors that share the ranges has been destroyed, then the program will crash accessing `ranges`.
1090: Fortran Note:
1091: The argument `ranges` must be declared as
1092: .vb
1093: PetscInt, pointer :: ranges(:)
1094: .ve
1095: and you have to return it with a call to `VecRestoreOwnershipRanges()` when no longer needed
1097: .seealso: [](ch_vectors), `Vec`, `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `VecGetOwnershipRange()`, `PetscSplitOwnership()`,
1098: `VecSetSizes()`, `VecCreateMPI()`, `PetscLayout`, `DMDAGetGhostCorners()`, `DM`
1099: @*/
1100: PetscErrorCode VecGetOwnershipRanges(Vec x, const PetscInt *ranges[])
1101: {
1102: PetscFunctionBegin;
1105: PetscCall(PetscLayoutGetRanges(x->map, ranges));
1106: PetscFunctionReturn(PETSC_SUCCESS);
1107: }
1109: /*@
1110: VecSetOption - Sets an option for controlling a vector's behavior with `VecSetValues()` and related routines
1112: Collective
1114: Input Parameters:
1115: + x - the vector
1116: . op - the `VecOption`
1117: - flag - turn the option on or off
1119: Level: intermediate
1121: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`, `VecOption`, `MatSetOption()`
1122: @*/
1123: PetscErrorCode VecSetOption(Vec x, VecOption op, PetscBool flag)
1124: {
1125: PetscFunctionBegin;
1128: PetscTryTypeMethod(x, setoption, op, flag);
1129: PetscFunctionReturn(PETSC_SUCCESS);
1130: }
1132: /* Default routines for obtaining and releasing; */
1133: /* may be used by any implementation */
1134: PetscErrorCode VecDuplicateVecs_Default(Vec w, PetscInt m, Vec *V[])
1135: {
1136: PetscFunctionBegin;
1137: PetscCheck(m > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "m must be > 0: m = %" PetscInt_FMT, m);
1138: PetscCall(PetscMalloc1(m, V));
1139: for (PetscInt i = 0; i < m; i++) PetscCall(VecDuplicate(w, *V + i));
1140: PetscFunctionReturn(PETSC_SUCCESS);
1141: }
1143: PetscErrorCode VecDestroyVecs_Default(PetscInt m, Vec v[])
1144: {
1145: PetscFunctionBegin;
1146: PetscAssertPointer(v, 2);
1147: for (PetscInt i = 0; i < m; i++) PetscCall(VecDestroy(&v[i]));
1148: PetscCall(PetscFree(v));
1149: PetscFunctionReturn(PETSC_SUCCESS);
1150: }
1152: /*@
1153: VecResetArray - Resets a vector to use its default memory. Call this
1154: after the use of `VecPlaceArray()`.
1156: Not Collective
1158: Input Parameter:
1159: . vec - the vector
1161: Level: developer
1163: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecPlaceArray()`
1164: @*/
1165: PetscErrorCode VecResetArray(Vec vec)
1166: {
1167: PetscFunctionBegin;
1170: PetscUseTypeMethod(vec, resetarray);
1171: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
1172: PetscFunctionReturn(PETSC_SUCCESS);
1173: }
1175: /*@
1176: VecLoad - Loads a vector that has been stored in binary or HDF5 format
1177: with `VecView()`.
1179: Collective
1181: Input Parameters:
1182: + vec - the newly loaded vector, this needs to have been created with `VecCreate()` or
1183: some related function before the call to `VecLoad()`.
1184: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()` or
1185: HDF5 file viewer, obtained from `PetscViewerHDF5Open()`
1187: Level: intermediate
1189: Notes:
1190: Defaults to the standard `VECSEQ` or `VECMPI`, if you want some other type of `Vec` call `VecSetFromOptions()`
1191: before calling this.
1193: The input file must contain the full global vector, as
1194: written by the routine `VecView()`.
1196: If the type or size of `vec` is not set before a call to `VecLoad()`, PETSc
1197: sets the type and the local and global sizes based on the vector it is reading in. If type and/or
1198: sizes are already set, then the same are used.
1200: If using the binary viewer and the blocksize of the vector is greater than one then you must provide a unique prefix to
1201: the vector with `PetscObjectSetOptionsPrefix`((`PetscObject`)vec,"uniqueprefix"); BEFORE calling `VecView()` on the
1202: vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
1203: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
1204: filename. If you copy the binary file, make sure you copy the associated .info file with it.
1206: If using HDF5, you must assign the `Vec` the same name as was used in the Vec
1207: that was stored in the file using `PetscObjectSetName()`. Otherwise you will
1208: get the error message: "Cannot H5DOpen2() with `Vec` name NAMEOFOBJECT".
1210: If the HDF5 file contains a two dimensional array the first dimension is treated as the block size
1211: in loading the vector. Hence, for example, using MATLAB notation h5create('vector.dat','/Test_Vec',[27 1]);
1212: will load a vector of size 27 and block size 27 thus resulting in all 27 entries being on the first process of
1213: vectors communicator and the rest of the processes having zero entries
1215: Notes for advanced users when using the binary viewer:
1216: Most users should not need to know the details of the binary storage
1217: format, since `VecLoad()` and `VecView()` completely hide these details.
1218: But for anyone who's interested, the standard binary vector storage
1219: format is
1220: .vb
1221: PetscInt VEC_FILE_CLASSID
1222: PetscInt number of rows
1223: PetscScalar *values of all entries
1224: .ve
1226: In addition, PETSc automatically uses byte swapping to work on all machines; the files
1227: are written ALWAYS using big-endian ordering. On small-endian machines the numbers
1228: are converted to the small-endian format when they are read in from the file.
1229: See PetscBinaryRead() and PetscBinaryWrite() to see how this may be done.
1231: .seealso: [](ch_vectors), `Vec`, `PetscViewerBinaryOpen()`, `VecView()`, `MatLoad()`
1232: @*/
1233: PetscErrorCode VecLoad(Vec vec, PetscViewer viewer)
1234: {
1235: PetscViewerFormat format;
1237: PetscFunctionBegin;
1240: PetscCheckSameComm(vec, 1, viewer, 2);
1242: PetscCall(VecSetErrorIfLocked(vec, 1));
1243: if (!((PetscObject)vec)->type_name && !vec->ops->create) PetscCall(VecSetType(vec, VECSTANDARD));
1244: PetscCall(PetscLogEventBegin(VEC_Load, viewer, 0, 0, 0));
1245: PetscCall(PetscViewerGetFormat(viewer, &format));
1246: if (format == PETSC_VIEWER_NATIVE && vec->ops->loadnative) {
1247: PetscUseTypeMethod(vec, loadnative, viewer);
1248: } else {
1249: PetscUseTypeMethod(vec, load, viewer);
1250: }
1251: PetscCall(PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0));
1252: PetscFunctionReturn(PETSC_SUCCESS);
1253: }
1255: /*@
1256: VecReciprocal - Replaces each component of a vector by its reciprocal.
1258: Logically Collective
1260: Input Parameter:
1261: . vec - the vector
1263: Output Parameter:
1264: . vec - the vector reciprocal
1266: Level: intermediate
1268: Note:
1269: Vector entries with value 0.0 are not changed
1271: .seealso: [](ch_vectors), `Vec`, `VecLog()`, `VecExp()`, `VecSqrtAbs()`
1272: @*/
1273: PetscErrorCode VecReciprocal(Vec vec)
1274: {
1275: PetscFunctionBegin;
1276: PetscCall(VecReciprocalAsync_Private(vec, NULL));
1277: PetscFunctionReturn(PETSC_SUCCESS);
1278: }
1280: /*@
1281: VecSetOperation - Allows the user to override a particular vector operation.
1283: Logically Collective; No Fortran Support
1285: Input Parameters:
1286: + vec - The vector to modify
1287: . op - The name of the operation
1288: - f - The function that provides the operation.
1290: Level: advanced
1292: Example Usage:
1293: .vb
1294: // some new VecView() implementation, must have the same signature as the function it seeks
1295: // to replace
1296: PetscErrorCode UserVecView(Vec x, PetscViewer viewer)
1297: {
1298: PetscFunctionBeginUser;
1299: // ...
1300: PetscFunctionReturn(PETSC_SUCCESS);
1301: }
1303: // Create a VECMPI which has a pre-defined VecView() implementation
1304: VecCreateMPI(comm, n, N, &x);
1305: // Calls the VECMPI implementation for VecView()
1306: VecView(x, viewer);
1308: VecSetOperation(x, VECOP_VIEW, (PetscErrorCodeFn *)UserVecView);
1309: // Now calls UserVecView()
1310: VecView(x, viewer);
1311: .ve
1313: Notes:
1314: `f` may be `NULL` to remove the operation from `vec`. Depending on the operation this may be
1315: allowed, however some always expect a valid function. In these cases an error will be raised
1316: when calling the interface routine in question.
1318: See `VecOperation` for an up-to-date list of override-able operations. The operations listed
1319: there have the form `VECOP_<OPERATION>`, where `<OPERATION>` is the suffix (in all capital
1320: letters) of the public interface routine (e.g., `VecView()` -> `VECOP_VIEW`).
1322: Overriding a particular `Vec`'s operation has no affect on any other `Vec`s past, present,
1323: or future. The user should also note that overriding a method is "destructive"; the previous
1324: method is not retained in any way.
1326: Each function MUST return `PETSC_SUCCESS` on success and
1327: nonzero on failure.
1329: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecGetOperation()`, `MatSetOperation()`, `MatShellSetOperation()`
1330: @*/
1331: PetscErrorCode VecSetOperation(Vec vec, VecOperation op, PetscErrorCodeFn *f)
1332: {
1333: PetscFunctionBegin;
1335: if (op == VECOP_VIEW && !vec->ops->viewnative) {
1336: vec->ops->viewnative = vec->ops->view;
1337: } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
1338: vec->ops->loadnative = vec->ops->load;
1339: }
1340: ((PetscErrorCodeFn **)vec->ops)[(int)op] = f;
1341: PetscFunctionReturn(PETSC_SUCCESS);
1342: }
1344: /*@
1345: VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1346: used during the assembly process to store values that belong to
1347: other processors.
1349: Not Collective, different processes can have different size stashes
1351: Input Parameters:
1352: + vec - the vector
1353: . size - the initial size of the stash.
1354: - bsize - the initial size of the block-stash(if used).
1356: Options Database Keys:
1357: + -vecstash_initial_size size or size0,size1,...,sizep-1 - set initial size
1358: - -vecstash_block_initial_size bsize or bsize0,bsize1,...,bsizep-1 - set initial block size
1360: Level: intermediate
1362: Notes:
1363: The block-stash is used for values set with `VecSetValuesBlocked()` while
1364: the stash is used for values set with `VecSetValues()`
1366: Run with the option -info and look for output of the form
1367: VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1368: to determine the appropriate value, MM, to use for size and
1369: VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1370: to determine the value, BMM to use for bsize
1372: PETSc attempts to smartly manage the stash size so there is little likelihood setting a
1373: a specific value here will affect performance
1375: .seealso: [](ch_vectors), `Vec`, `VecSetBlockSize()`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecStashView()`
1376: @*/
1377: PetscErrorCode VecStashSetInitialSize(Vec vec, PetscInt size, PetscInt bsize)
1378: {
1379: PetscFunctionBegin;
1381: PetscCall(VecStashSetInitialSize_Private(&vec->stash, size));
1382: PetscCall(VecStashSetInitialSize_Private(&vec->bstash, bsize));
1383: PetscFunctionReturn(PETSC_SUCCESS);
1384: }
1386: /*@
1387: VecSetRandom - Sets all components of a vector to random numbers.
1389: Logically Collective
1391: Input Parameters:
1392: + x - the vector
1393: - rctx - the random number context, formed by `PetscRandomCreate()`, or use `NULL` and it will create one internally.
1395: Output Parameter:
1396: . x - the vector
1398: Example of Usage:
1399: .vb
1400: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1401: VecSetRandom(x,rctx);
1402: PetscRandomDestroy(&rctx);
1403: .ve
1405: Level: intermediate
1407: .seealso: [](ch_vectors), `Vec`, `VecSet()`, `VecSetValues()`, `PetscRandomCreate()`, `PetscRandomDestroy()`
1408: @*/
1409: PetscErrorCode VecSetRandom(Vec x, PetscRandom rctx)
1410: {
1411: PetscRandom randObj = NULL;
1413: PetscFunctionBegin;
1417: VecCheckAssembled(x);
1418: PetscCall(VecSetErrorIfLocked(x, 1));
1420: if (!rctx) {
1421: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)x), &randObj));
1422: PetscCall(PetscRandomSetType(randObj, x->defaultrandtype));
1423: PetscCall(PetscRandomSetFromOptions(randObj));
1424: rctx = randObj;
1425: }
1427: PetscCall(PetscLogEventBegin(VEC_SetRandom, x, rctx, 0, 0));
1428: PetscUseTypeMethod(x, setrandom, rctx);
1429: PetscCall(PetscLogEventEnd(VEC_SetRandom, x, rctx, 0, 0));
1431: PetscCall(PetscRandomDestroy(&randObj));
1432: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1433: PetscFunctionReturn(PETSC_SUCCESS);
1434: }
1436: /*@
1437: VecSetRandomGaussian - Fills a vector with Gaussian random values of the given mean and standard deviation.
1439: Collective
1441: Input Parameters:
1442: + v - the vector to fill
1443: . rng - PETSc random number generator
1444: . mean - desired mean of the Gaussian samples
1445: - std_dev - desired standard deviation
1447: Level: advanced
1449: Note:
1450: For complex builds where `PetscScalar` is complex the imaginary part of all the vector entries is zero
1452: Developer Note:
1453: Uses the Box-Muller transform to generate normally distributed random numbers
1454: from uniform random numbers. Handles edge cases where uniform random values
1455: approach 0 or 1.
1457: .seealso: [](ch_vectors), [](ch_da), `PetscDA`, `PetscRandom`, `PetscRandomSetInterval()`, `VecSetRandom()`
1458: @*/
1459: PetscErrorCode VecSetRandomGaussian(Vec v, PetscRandom rng, PetscReal mean, PetscReal std_dev)
1460: {
1461: PetscInt n;
1462: PetscScalar *array;
1463: PetscReal u1, u2;
1464: PetscReal gauss_sample1, gauss_sample2, magnitude, theta;
1465: const PetscReal min_uniform = PETSC_MACHINE_EPSILON;
1466: const PetscInt max_retry_count = 100;
1468: PetscFunctionBegin;
1473: PetscCheck(!PetscIsInfOrNanReal(mean), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Mean must be a finite real number");
1474: PetscCheck(std_dev >= 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Standard deviation must be non-negative, got %g", (double)std_dev);
1475: PetscCheck(!PetscIsInfOrNanReal(std_dev), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Standard deviation must be a finite real number");
1477: PetscCall(VecGetLocalSize(v, &n));
1478: if (n == 0) PetscFunctionReturn(PETSC_SUCCESS);
1480: if (std_dev == 0.0) {
1481: PetscCall(VecSet(v, mean));
1482: PetscFunctionReturn(PETSC_SUCCESS);
1483: }
1485: PetscCall(VecGetArrayWrite(v, &array));
1487: /*
1488: Generate Gaussian-distributed random values using the Box-Muller transform.
1489: This transform converts pairs of uniform random variables U1, U2 ~ Uniform(0,1)
1490: into pairs of independent standard normal variables Z0, Z1 ~ N(0,1):
1491: Z0 = sqrt(-2 * ln(U1)) * cos(2pi * U2)
1492: Z1 = sqrt(-2 * ln(U1)) * sin(2pi * U2)
1493: Then scale and shift to get desired mean and standard deviation.
1494: */
1495: for (PetscInt i = 0; i < n; i += 2) {
1496: PetscInt retry_count = 0;
1498: /*
1499: Generate U1 and ensure it's not too close to 0 to avoid log(0) singularity.
1500: Add retry limit to prevent infinite loops in case of RNG failure.
1501: */
1502: do {
1503: PetscCall(PetscRandomGetValueReal(rng, &u1));
1504: retry_count++;
1505: 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);
1506: } while (u1 < min_uniform);
1508: PetscCall(PetscRandomGetValueReal(rng, &u2));
1510: /*
1511: Apply Box-Muller transform:
1512: - magnitude: sqrt(-2 * ln(U1)) represents the radial distance from origin
1513: - theta: 2pi * U2 represents the angle uniformly distributed on [0, 2pi]
1514: - Converting from polar to Cartesian coordinates yields two independent samples
1515: */
1516: magnitude = PetscSqrtReal(-2.0 * PetscLogReal(u1));
1517: theta = 2.0 * PETSC_PI * u2;
1518: gauss_sample1 = magnitude * PetscCosReal(theta);
1519: gauss_sample2 = magnitude * PetscSinReal(theta);
1521: /* Scale and shift to achieve desired mean and standard deviation */
1522: array[i] = mean + std_dev * gauss_sample1;
1523: if (i + 1 < n) array[i + 1] = mean + std_dev * gauss_sample2;
1524: }
1526: PetscCall(VecRestoreArrayWrite(v, &array));
1527: PetscFunctionReturn(PETSC_SUCCESS);
1528: }
1530: /*@
1531: VecZeroEntries - puts a `0.0` in each element of a vector
1533: Logically Collective
1535: Input Parameter:
1536: . vec - The vector
1538: Level: beginner
1540: Note:
1541: If the norm of the vector is known to be zero then this skips the unneeded zeroing process
1543: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecSetOptionsPrefix()`, `VecSet()`, `VecSetValues()`
1544: @*/
1545: PetscErrorCode VecZeroEntries(Vec vec)
1546: {
1547: PetscFunctionBegin;
1548: PetscCall(VecSet(vec, 0));
1549: PetscFunctionReturn(PETSC_SUCCESS);
1550: }
1552: /*
1553: VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
1554: processor and a PETSc MPI vector on more than one processor.
1556: Collective
1558: Input Parameter:
1559: . vec - The vector
1561: Level: intermediate
1563: .seealso: [](ch_vectors), `Vec`, `VecSetFromOptions()`, `VecSetType()`
1564: */
1565: static PetscErrorCode VecSetTypeFromOptions_Private(Vec vec, PetscOptionItems PetscOptionsObject)
1566: {
1567: PetscBool opt;
1568: VecType defaultType;
1569: char typeName[256];
1570: PetscMPIInt size;
1572: PetscFunctionBegin;
1573: if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1574: else {
1575: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
1576: if (size > 1) defaultType = VECMPI;
1577: else defaultType = VECSEQ;
1578: }
1580: PetscCall(VecRegisterAll());
1581: PetscCall(PetscOptionsFList("-vec_type", "Vector type", "VecSetType", VecList, defaultType, typeName, sizeof(typeName), &opt));
1582: if (opt) PetscCall(VecSetType(vec, typeName));
1583: else PetscCall(VecSetType(vec, defaultType));
1584: PetscFunctionReturn(PETSC_SUCCESS);
1585: }
1587: /*@
1588: VecSetFromOptions - Configures the vector from the options database.
1590: Collective
1592: Input Parameter:
1593: . vec - The vector
1595: Options Database Key:
1596: . -vec_type type - set the vector type, see `VecType`
1598: Level: beginner
1600: Notes:
1601: To see all options, run your program with the -help option.
1603: Must be called after `VecCreate()` but before the vector is used.
1605: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecSetOptionsPrefix()`, `VecType`
1606: @*/
1607: PetscErrorCode VecSetFromOptions(Vec vec)
1608: {
1609: PetscBool flg;
1610: PetscInt bind_below = 0;
1612: PetscFunctionBegin;
1615: PetscObjectOptionsBegin((PetscObject)vec);
1616: /* Handle vector type options */
1617: PetscCall(VecSetTypeFromOptions_Private(vec, PetscOptionsObject));
1619: /* Handle specific vector options */
1620: PetscTryTypeMethod(vec, setfromoptions, PetscOptionsObject);
1622: /* Bind to CPU if below a user-specified size threshold.
1623: * This perhaps belongs in the options for the GPU Vec types, but VecBindToCPU() does nothing when called on non-GPU types,
1624: * and putting it here makes is more maintainable than duplicating this for all. */
1625: 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));
1626: if (flg && vec->map->n < bind_below) PetscCall(VecBindToCPU(vec, PETSC_TRUE));
1628: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1629: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)vec, PetscOptionsObject));
1630: PetscOptionsEnd();
1631: PetscFunctionReturn(PETSC_SUCCESS);
1632: }
1634: /*@
1635: VecSetSizes - Sets the local and global sizes, and checks to determine compatibility of the sizes
1637: Collective
1639: Input Parameters:
1640: + v - the vector
1641: . n - the local size (or `PETSC_DECIDE` to have it set)
1642: - N - the global size (or `PETSC_DETERMINE` to have it set)
1644: Level: intermediate
1646: Notes:
1647: `N` cannot be `PETSC_DETERMINE` if `n` is `PETSC_DECIDE`
1649: If one processor calls this with `N` of `PETSC_DETERMINE` then all processors must, otherwise the program will hang.
1651: If `n` is not `PETSC_DECIDE`, then the value determines the `PetscLayout` of the vector and the ranges returned by
1652: `VecGetOwnershipRange()` and `VecGetOwnershipRanges()`
1654: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecCreateSeq()`, `VecCreateMPI()`, `VecGetSize()`, `PetscSplitOwnership()`, `PetscLayout`,
1655: `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`, `MatSetSizes()`
1656: @*/
1657: PetscErrorCode VecSetSizes(Vec v, PetscInt n, PetscInt N)
1658: {
1659: PetscFunctionBegin;
1661: if (N >= 0) {
1663: PetscCheck(n <= N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local size %" PetscInt_FMT " cannot be larger than global size %" PetscInt_FMT, n, N);
1664: }
1665: 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,
1666: v->map->n, v->map->N);
1667: v->map->n = n;
1668: v->map->N = N;
1669: PetscTryTypeMethod(v, create);
1670: v->ops->create = NULL;
1671: PetscFunctionReturn(PETSC_SUCCESS);
1672: }
1674: /*@
1675: VecSetBlockSize - Sets the block size for future calls to `VecSetValuesBlocked()`
1676: and `VecSetValuesBlockedLocal()`.
1678: Logically Collective
1680: Input Parameters:
1681: + v - the vector
1682: - bs - the blocksize
1684: Level: advanced
1686: Note:
1687: All vectors obtained by `VecDuplicate()` inherit the same blocksize.
1689: Vectors obtained with `DMCreateGlobalVector()` and `DMCreateLocalVector()` generally already have a blocksize set based on the state of the `DM`
1691: .seealso: [](ch_vectors), `Vec`, `VecSetValuesBlocked()`, `VecSetLocalToGlobalMapping()`, `VecGetBlockSize()`
1692: @*/
1693: PetscErrorCode VecSetBlockSize(Vec v, PetscInt bs)
1694: {
1695: PetscFunctionBegin;
1698: PetscCall(PetscLayoutSetBlockSize(v->map, bs));
1699: v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1700: PetscFunctionReturn(PETSC_SUCCESS);
1701: }
1703: /*@
1704: VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for `VecSetValuesBlocked()`
1705: and `VecSetValuesBlockedLocal()`.
1707: Not Collective
1709: Input Parameter:
1710: . v - the vector
1712: Output Parameter:
1713: . bs - the blocksize
1715: Level: advanced
1717: Note:
1718: All vectors obtained by `VecDuplicate()` inherit the same blocksize.
1720: .seealso: [](ch_vectors), `Vec`, `VecSetValuesBlocked()`, `VecSetLocalToGlobalMapping()`, `VecSetBlockSize()`
1721: @*/
1722: PetscErrorCode VecGetBlockSize(Vec v, PetscInt *bs)
1723: {
1724: PetscFunctionBegin;
1726: PetscAssertPointer(bs, 2);
1727: PetscCall(PetscLayoutGetBlockSize(v->map, bs));
1728: PetscFunctionReturn(PETSC_SUCCESS);
1729: }
1731: /*@
1732: VecSetOptionsPrefix - Sets the prefix used for searching for all
1733: `Vec` options in the database.
1735: Logically Collective
1737: Input Parameters:
1738: + v - the `Vec` context
1739: - prefix - the prefix to prepend to all option names
1741: Level: advanced
1743: Note:
1744: A hyphen (-) must NOT be given at the beginning of the prefix name.
1745: The first character of all runtime options is AUTOMATICALLY the hyphen.
1747: .seealso: [](ch_vectors), `Vec`, `VecSetFromOptions()`
1748: @*/
1749: PetscErrorCode VecSetOptionsPrefix(Vec v, const char prefix[])
1750: {
1751: PetscFunctionBegin;
1753: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)v, prefix));
1754: PetscFunctionReturn(PETSC_SUCCESS);
1755: }
1757: /*@
1758: VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1759: `Vec` options in the database.
1761: Logically Collective
1763: Input Parameters:
1764: + v - the `Vec` context
1765: - prefix - the prefix to prepend to all option names
1767: Level: advanced
1769: Note:
1770: A hyphen (-) must NOT be given at the beginning of the prefix name.
1771: The first character of all runtime options is AUTOMATICALLY the hyphen.
1773: .seealso: [](ch_vectors), `Vec`, `VecGetOptionsPrefix()`
1774: @*/
1775: PetscErrorCode VecAppendOptionsPrefix(Vec v, const char prefix[])
1776: {
1777: PetscFunctionBegin;
1779: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)v, prefix));
1780: PetscFunctionReturn(PETSC_SUCCESS);
1781: }
1783: /*@
1784: VecGetOptionsPrefix - Sets the prefix used for searching for all
1785: Vec options in the database.
1787: Not Collective
1789: Input Parameter:
1790: . v - the `Vec` context
1792: Output Parameter:
1793: . prefix - pointer to the prefix string used
1795: Level: advanced
1797: .seealso: [](ch_vectors), `Vec`, `VecAppendOptionsPrefix()`
1798: @*/
1799: PetscErrorCode VecGetOptionsPrefix(Vec v, const char *prefix[])
1800: {
1801: PetscFunctionBegin;
1803: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)v, prefix));
1804: PetscFunctionReturn(PETSC_SUCCESS);
1805: }
1807: /*@
1808: VecGetState - Gets the state of a `Vec`.
1810: Not Collective
1812: Input Parameter:
1813: . v - the `Vec` context
1815: Output Parameter:
1816: . state - the object state
1818: Level: advanced
1820: Note:
1821: Object state is an integer which gets increased every time
1822: the object is changed. By saving and later querying the object state
1823: one can determine whether information about the object is still current.
1825: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `PetscObjectStateGet()`
1826: @*/
1827: PetscErrorCode VecGetState(Vec v, PetscObjectState *state)
1828: {
1829: PetscFunctionBegin;
1831: PetscAssertPointer(state, 2);
1832: PetscCall(PetscObjectStateGet((PetscObject)v, state));
1833: PetscFunctionReturn(PETSC_SUCCESS);
1834: }
1836: /*@
1837: VecSetUp - Sets up the internal vector data structures for the later use.
1839: Collective
1841: Input Parameter:
1842: . v - the `Vec` context
1844: Level: advanced
1846: Notes:
1847: For basic use of the `Vec` classes the user need not explicitly call
1848: `VecSetUp()`, since these actions will happen automatically.
1850: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecDestroy()`
1851: @*/
1852: PetscErrorCode VecSetUp(Vec v)
1853: {
1854: PetscMPIInt size;
1856: PetscFunctionBegin;
1858: PetscCheck(v->map->n >= 0 || v->map->N >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Sizes not set");
1859: if (!((PetscObject)v)->type_name) {
1860: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1861: if (size == 1) PetscCall(VecSetType(v, VECSEQ));
1862: else PetscCall(VecSetType(v, VECMPI));
1863: }
1864: PetscFunctionReturn(PETSC_SUCCESS);
1865: }
1867: /*
1868: These currently expose the PetscScalar/PetscReal in updating the
1869: cached norm. If we push those down into the implementation these
1870: will become independent of PetscScalar/PetscReal
1871: */
1873: PetscErrorCode VecCopyAsync_Private(Vec x, Vec y, PetscDeviceContext dctx)
1874: {
1875: PetscBool flgs[4];
1876: PetscReal norms[4] = {0.0, 0.0, 0.0, 0.0};
1878: PetscFunctionBegin;
1883: if (x == y) PetscFunctionReturn(PETSC_SUCCESS);
1884: VecCheckSameLocalSize(x, 1, y, 2);
1885: VecCheckAssembled(x);
1886: PetscCall(VecSetErrorIfLocked(y, 2));
1888: #if !PetscDefined(USE_MIXED_PRECISION)
1889: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
1890: #endif
1892: PetscCall(PetscLogEventBegin(VEC_Copy, x, y, 0, 0));
1893: #if PetscDefined(USE_MIXED_PRECISION)
1894: extern PetscErrorCode VecGetArray(Vec, double **);
1895: extern PetscErrorCode VecRestoreArray(Vec, double **);
1896: extern PetscErrorCode VecGetArray(Vec, float **);
1897: extern PetscErrorCode VecRestoreArray(Vec, float **);
1898: extern PetscErrorCode VecGetArrayRead(Vec, const double **);
1899: extern PetscErrorCode VecRestoreArrayRead(Vec, const double **);
1900: extern PetscErrorCode VecGetArrayRead(Vec, const float **);
1901: extern PetscErrorCode VecRestoreArrayRead(Vec, const float **);
1902: if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1903: PetscInt i, n;
1904: const float *xx;
1905: double *yy;
1906: PetscCall(VecGetArrayRead(x, &xx));
1907: PetscCall(VecGetArray(y, &yy));
1908: PetscCall(VecGetLocalSize(x, &n));
1909: for (i = 0; i < n; i++) yy[i] = xx[i];
1910: PetscCall(VecRestoreArrayRead(x, &xx));
1911: PetscCall(VecRestoreArray(y, &yy));
1912: } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1913: PetscInt i, n;
1914: float *yy;
1915: const double *xx;
1916: PetscCall(VecGetArrayRead(x, &xx));
1917: PetscCall(VecGetArray(y, &yy));
1918: PetscCall(VecGetLocalSize(x, &n));
1919: for (i = 0; i < n; i++) yy[i] = (float)xx[i];
1920: PetscCall(VecRestoreArrayRead(x, &xx));
1921: PetscCall(VecRestoreArray(y, &yy));
1922: } else PetscUseTypeMethod(x, copy, y);
1923: #else
1924: VecMethodDispatch(x, dctx, VecAsyncFnName(Copy), copy, (Vec, Vec, PetscDeviceContext), y);
1925: #endif
1927: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1928: #if !PetscDefined(USE_MIXED_PRECISION)
1929: for (PetscInt i = 0; i < 4; i++) {
1930: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)y, NormIds[i], norms[i]));
1931: }
1932: #endif
1934: PetscCall(PetscLogEventEnd(VEC_Copy, x, y, 0, 0));
1935: PetscFunctionReturn(PETSC_SUCCESS);
1936: }
1938: /*@
1939: VecCopy - Copies a vector `y = x`
1941: Logically Collective
1943: Input Parameter:
1944: . x - the vector
1946: Output Parameter:
1947: . y - the copy
1949: Level: beginner
1951: Note:
1952: For default parallel PETSc vectors, both `x` and `y` must be distributed in
1953: the same manner; local copies are done.
1955: Developer Notes:
1956: `PetscCheckSameTypeAndComm`(x,1,y,2) is not used on these vectors because we allow one
1957: of the vectors to be sequential and one to be parallel so long as both have the same
1958: local sizes. This is used in some internal functions in PETSc.
1960: .seealso: [](ch_vectors), `Vec`, `VecDuplicate()`
1961: @*/
1962: PetscErrorCode VecCopy(Vec x, Vec y)
1963: {
1964: PetscFunctionBegin;
1965: PetscCall(VecCopyAsync_Private(x, y, NULL));
1966: PetscFunctionReturn(PETSC_SUCCESS);
1967: }
1969: PetscErrorCode VecSwapAsync_Private(Vec x, Vec y, PetscDeviceContext dctx)
1970: {
1971: PetscReal normxs[4], normys[4];
1972: PetscBool flgxs[4], flgys[4];
1974: PetscFunctionBegin;
1979: PetscCheckSameTypeAndComm(x, 1, y, 2);
1980: VecCheckSameSize(x, 1, y, 2);
1981: VecCheckAssembled(x);
1982: VecCheckAssembled(y);
1983: PetscCall(VecSetErrorIfLocked(x, 1));
1984: PetscCall(VecSetErrorIfLocked(y, 2));
1986: for (PetscInt i = 0; i < 4; i++) {
1987: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], normxs[i], flgxs[i]));
1988: PetscCall(PetscObjectComposedDataGetReal((PetscObject)y, NormIds[i], normys[i], flgys[i]));
1989: }
1991: PetscCall(PetscLogEventBegin(VEC_Swap, x, y, 0, 0));
1992: VecMethodDispatch(x, dctx, VecAsyncFnName(Swap), swap, (Vec, Vec, PetscDeviceContext), y);
1993: PetscCall(PetscLogEventEnd(VEC_Swap, x, y, 0, 0));
1995: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1996: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1997: for (PetscInt i = 0; i < 4; i++) {
1998: if (flgxs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)y, NormIds[i], normxs[i]));
1999: if (flgys[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], normys[i]));
2000: }
2001: PetscFunctionReturn(PETSC_SUCCESS);
2002: }
2003: /*@
2004: VecSwap - Swaps the values between two vectors, `x` and `y`.
2006: Logically Collective
2008: Input Parameters:
2009: + x - the first vector
2010: - y - the second vector
2012: Level: advanced
2014: .seealso: [](ch_vectors), `Vec`, `VecSet()`
2015: @*/
2016: PetscErrorCode VecSwap(Vec x, Vec y)
2017: {
2018: PetscFunctionBegin;
2019: PetscCall(VecSwapAsync_Private(x, y, NULL));
2020: PetscFunctionReturn(PETSC_SUCCESS);
2021: }
2023: /*@
2024: VecStashViewFromOptions - Processes command line options to determine if/how a `VecStash` object is to be viewed.
2026: Collective
2028: Input Parameters:
2029: + obj - the `Vec` containing a stash
2030: . bobj - optional other object that provides the options prefix, pass `NULL` to use the options prefix of `obj`
2031: - name - option to activate viewing
2033: Options Database Key:
2034: . -name viewer_specification - See `PetscOptionsCreateViewer()` for the values of `viewer_specification`
2036: Level: intermediate
2038: Note:
2039: This checks the options database, creates the viewer on-the-fly, uses it and then destroys it. Hence it should not be called in heavily used routines,
2040: rather `PetscOptionsCreateViewer()` should be used to construct the viewer once which can then be utilized in the heavily used routine.
2042: Developer Notes:
2043: This cannot use `PetscObjectViewFromOptions()` because it takes a `Vec` as an argument but does not use `VecView()`
2045: .seealso: [](ch_vectors), `Vec`, `VecStashView()`, `VecStashSetInitialSize()`, `PetscOptionsCreateViewer()`
2046: @*/
2047: PetscErrorCode VecStashViewFromOptions(Vec obj, PetscObject bobj, const char name[])
2048: {
2049: PetscViewer viewer;
2050: PetscBool flg;
2051: PetscViewerFormat format;
2052: char *prefix;
2054: PetscFunctionBegin;
2055: prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
2056: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)obj), ((PetscObject)obj)->options, prefix, name, &viewer, &format, &flg));
2057: if (flg) {
2058: PetscCall(PetscViewerPushFormat(viewer, format));
2059: PetscCall(VecStashView(obj, viewer));
2060: PetscCall(PetscViewerPopFormat(viewer));
2061: PetscCall(PetscViewerDestroy(&viewer));
2062: }
2063: PetscFunctionReturn(PETSC_SUCCESS);
2064: }
2066: /*@
2067: VecStashView - Prints the entries in the vector stash and block stash.
2069: Collective
2071: Input Parameters:
2072: + v - the vector
2073: - viewer - the viewer
2075: Level: advanced
2077: .seealso: [](ch_vectors), `Vec`, `VecSetBlockSize()`, `VecSetValues()`, `VecSetValuesBlocked()`
2078: @*/
2079: PetscErrorCode VecStashView(Vec v, PetscViewer viewer)
2080: {
2081: PetscMPIInt rank;
2082: PetscInt i;
2083: PetscBool match;
2084: VecStash *s;
2085: PetscScalar val;
2087: PetscFunctionBegin;
2090: PetscCheckSameComm(v, 1, viewer, 2);
2092: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &match));
2093: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_SUP, "Stash viewer only works with ASCII viewer not %s", ((PetscObject)v)->type_name);
2094: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
2095: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank));
2096: s = &v->bstash;
2098: /* print block stash */
2099: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2100: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Vector Block stash size %" PetscInt_FMT " block size %" PetscInt_FMT "\n", rank, s->n, s->bs));
2101: for (i = 0; i < s->n; i++) {
2102: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " ", rank, s->idx[i]));
2103: for (PetscInt j = 0; j < s->bs; j++) {
2104: val = s->array[i * s->bs + j];
2105: #if PetscDefined(USE_COMPLEX)
2106: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "(%18.16e %18.16e) ", (double)PetscRealPart(val), (double)PetscImaginaryPart(val)));
2107: #else
2108: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%18.16e ", (double)val));
2109: #endif
2110: }
2111: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2112: }
2113: PetscCall(PetscViewerFlush(viewer));
2115: s = &v->stash;
2117: /* print basic stash */
2118: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Vector stash size %" PetscInt_FMT "\n", rank, s->n));
2119: for (i = 0; i < s->n; i++) {
2120: val = s->array[i];
2121: #if PetscDefined(USE_COMPLEX)
2122: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " (%18.16e %18.16e) ", rank, s->idx[i], (double)PetscRealPart(val), (double)PetscImaginaryPart(val)));
2123: #else
2124: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " %18.16e\n", rank, s->idx[i], (double)val));
2125: #endif
2126: }
2127: PetscCall(PetscViewerFlush(viewer));
2128: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2129: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
2130: PetscFunctionReturn(PETSC_SUCCESS);
2131: }
2133: /*@
2134: PetscOptionsGetVec - Gets a `Vec` from the options database as an array of real values
2136: Collective
2138: Input Parameters:
2139: + options - the options database, or `NULL` for the default global one
2140: . prefix - an option prefix, or `NULL`
2141: - key - the option name (must include the leading `-`)
2143: Output Parameters:
2144: + v - the vector to fill in on option match; unchanged if the option is not found
2145: - set - `PETSC_TRUE` if the option was found (may be `NULL`)
2147: Level: intermediate
2149: Note:
2150: The option value is read as an array of `PetscReal` of length equal to the global size of `v`; each MPI process
2151: copies the entries corresponding to its local ownership range into `v`.
2153: .seealso: `Vec`, `PetscOptionsGetRealArray()`, `PetscOptionsGetInt()`, `PetscOptionsGetReal()`, `VecView()`
2154: @*/
2155: PetscErrorCode PetscOptionsGetVec(PetscOptions options, const char prefix[], const char key[], Vec v, PetscBool *set)
2156: {
2157: PetscInt i, N, rstart, rend;
2158: PetscScalar *xx;
2159: PetscReal *xreal;
2160: PetscBool iset;
2162: PetscFunctionBegin;
2163: PetscCall(VecGetOwnershipRange(v, &rstart, &rend));
2164: PetscCall(VecGetSize(v, &N));
2165: PetscCall(PetscCalloc1(N, &xreal));
2166: PetscCall(PetscOptionsGetRealArray(options, prefix, key, xreal, &N, &iset));
2167: if (iset) {
2168: PetscCall(VecGetArray(v, &xx));
2169: for (i = rstart; i < rend; i++) xx[i - rstart] = xreal[i];
2170: PetscCall(VecRestoreArray(v, &xx));
2171: }
2172: PetscCall(PetscFree(xreal));
2173: if (set) *set = iset;
2174: PetscFunctionReturn(PETSC_SUCCESS);
2175: }
2177: /*@
2178: VecGetLayout - get `PetscLayout` describing a vector layout
2180: Not Collective
2182: Input Parameter:
2183: . x - the vector
2185: Output Parameter:
2186: . map - the layout
2188: Level: developer
2190: Note:
2191: The layout determines what vector elements are contained on each MPI process
2193: .seealso: [](ch_vectors), `PetscLayout`, `Vec`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
2194: @*/
2195: PetscErrorCode VecGetLayout(Vec x, PetscLayout *map)
2196: {
2197: PetscFunctionBegin;
2199: PetscAssertPointer(map, 2);
2200: *map = x->map;
2201: PetscFunctionReturn(PETSC_SUCCESS);
2202: }
2204: /*@
2205: VecSetLayout - set `PetscLayout` describing vector layout
2207: Not Collective
2209: Input Parameters:
2210: + x - the vector
2211: - map - the layout
2213: Level: developer
2215: Note:
2216: It is normally only valid to replace the layout with a layout known to be equivalent.
2218: .seealso: [](ch_vectors), `Vec`, `PetscLayout`, `VecGetLayout()`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
2219: @*/
2220: PetscErrorCode VecSetLayout(Vec x, PetscLayout map)
2221: {
2222: PetscFunctionBegin;
2224: PetscCall(PetscLayoutReference(map, &x->map));
2225: PetscFunctionReturn(PETSC_SUCCESS);
2226: }
2228: /*@
2229: VecFlag - set infinity into the local part of the vector on any subset of MPI processes
2231: Logically Collective
2233: Input Parameters:
2234: + xin - the vector, can be `NULL` but only if on all processes
2235: - flg - indicates if this processes portion of the vector should be set to infinity
2237: Level: developer
2239: Note:
2240: This removes the values from the vector norm cache for all processes by calling `PetscObjectIncrease()`.
2242: This is used for any subset of MPI processes to indicate an failure in a solver, after the next use of `VecNorm()` if
2243: `KSPCheckNorm()` detects an infinity and at least one of the MPI processes has a not converged reason then the `KSP`
2244: object collectively is labeled as not converged.
2246: .seealso: [](ch_vectors), `Vec`, `PetscLayout`, `VecGetLayout()`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`, `MatSetInf()`
2247: @*/
2248: PetscErrorCode VecFlag(Vec xin, PetscInt flg)
2249: {
2250: // MSVC gives "divide by zero" error at compile time - so declare as volatile to skip this check.
2251: volatile PetscReal one = 1.0, zero = 0.0;
2252: PetscScalar inf;
2254: PetscFunctionBegin;
2255: if (!xin) PetscFunctionReturn(PETSC_SUCCESS);
2257: PetscCall(PetscObjectStateIncrease((PetscObject)xin));
2258: if (flg) {
2259: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2260: inf = one / zero;
2261: PetscCall(PetscFPTrapPop());
2262: if (xin->ops->set) PetscUseTypeMethod(xin, set, inf);
2263: else {
2264: PetscInt n;
2265: PetscScalar *xx;
2267: PetscCall(VecGetLocalSize(xin, &n));
2268: PetscCall(VecGetArrayWrite(xin, &xx));
2269: for (PetscInt i = 0; i < n; ++i) xx[i] = inf;
2270: PetscCall(VecRestoreArrayWrite(xin, &xx));
2271: }
2272: }
2273: PetscFunctionReturn(PETSC_SUCCESS);
2274: }
2276: /*@
2277: VecSetInf - set infinity into the local part of the vector
2279: Not Collective
2281: Input Parameters:
2282: . xin - the vector
2284: Level: developer
2286: Note:
2287: Deprecated, see `VecFlag()`
2288: This is used for any subset of MPI processes to indicate an failure in a solver, after the next use of `VecNorm()` if
2289: `KSPCheckNorm()` detects an infinity and at least one of the MPI processes has a not converged reason then the `KSP`
2290: object collectively is labeled as not converged.
2292: This cannot be called if `xin` has a cached norm available
2294: .seealso: [](ch_vectors), `VecFlag()`, `Vec`, `PetscLayout`, `VecGetLayout()`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
2295: @*/
2296: PetscErrorCode VecSetInf(Vec xin)
2297: {
2298: // MSVC gives "divide by zero" error at compile time - so declare as volatile to skip this check.
2299: volatile PetscReal one = 1.0, zero = 0.0;
2300: PetscScalar inf;
2301: PetscBool flg;
2303: PetscFunctionBegin;
2304: PetscCall(VecNormAvailable(xin, NORM_2, &flg, NULL));
2305: PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot call VecSetInf() if the vector has a cached norm");
2306: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2307: inf = one / zero;
2308: PetscCall(PetscFPTrapPop());
2309: if (xin->ops->set) PetscUseTypeMethod(xin, set, inf);
2310: else {
2311: PetscInt n;
2312: PetscScalar *xx;
2314: PetscCall(VecGetLocalSize(xin, &n));
2315: PetscCall(VecGetArrayWrite(xin, &xx));
2316: for (PetscInt i = 0; i < n; ++i) xx[i] = inf;
2317: PetscCall(VecRestoreArrayWrite(xin, &xx));
2318: }
2319: PetscFunctionReturn(PETSC_SUCCESS);
2320: }
2322: /*@
2323: VecBindToCPU - marks a vector to temporarily stay on the CPU and perform computations on the CPU
2325: Logically collective
2327: Input Parameters:
2328: + v - the vector
2329: - flg - bind to the CPU if value of `PETSC_TRUE`
2331: Level: intermediate
2333: .seealso: [](ch_vectors), `Vec`, `VecBoundToCPU()`
2334: @*/
2335: PetscErrorCode VecBindToCPU(Vec v, PetscBool flg)
2336: {
2337: PetscFunctionBegin;
2340: #if PetscDefined(HAVE_DEVICE)
2341: if (v->boundtocpu == flg) PetscFunctionReturn(PETSC_SUCCESS);
2342: v->boundtocpu = flg;
2343: PetscTryTypeMethod(v, bindtocpu, flg);
2344: #endif
2345: PetscFunctionReturn(PETSC_SUCCESS);
2346: }
2348: /*@
2349: VecBoundToCPU - query if a vector is bound to the CPU
2351: Not collective
2353: Input Parameter:
2354: . v - the vector
2356: Output Parameter:
2357: . flg - the logical flag
2359: Level: intermediate
2361: .seealso: [](ch_vectors), `Vec`, `VecBindToCPU()`
2362: @*/
2363: PetscErrorCode VecBoundToCPU(Vec v, PetscBool *flg)
2364: {
2365: PetscFunctionBegin;
2367: PetscAssertPointer(flg, 2);
2368: #if PetscDefined(HAVE_DEVICE)
2369: *flg = v->boundtocpu;
2370: #else
2371: *flg = PETSC_TRUE;
2372: #endif
2373: PetscFunctionReturn(PETSC_SUCCESS);
2374: }
2376: /*@
2377: VecSetBindingPropagates - Sets whether the state of being bound to the CPU for a GPU vector type propagates to child and some other associated objects
2379: Input Parameters:
2380: + v - the vector
2381: - flg - flag indicating whether the boundtocpu flag should be propagated
2383: Level: developer
2385: Notes:
2386: 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.
2387: The created vectors will also have their bindingpropagates flag set to true.
2389: Developer Notes:
2390: If a `DMDA` has the `-dm_bind_below option` set to true, then vectors created by `DMCreateGlobalVector()` will have `VecSetBindingPropagates()` called on them to
2391: set their bindingpropagates flag to true.
2393: .seealso: [](ch_vectors), `Vec`, `MatSetBindingPropagates()`, `VecGetBindingPropagates()`
2394: @*/
2395: PetscErrorCode VecSetBindingPropagates(Vec v, PetscBool flg)
2396: {
2397: PetscFunctionBegin;
2399: #if PetscDefined(HAVE_VIENNACL) || PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
2400: v->bindingpropagates = flg;
2401: #endif
2402: PetscFunctionReturn(PETSC_SUCCESS);
2403: }
2405: /*@
2406: VecGetBindingPropagates - Gets whether the state of being bound to the CPU for a GPU vector type propagates to child and some other associated objects
2408: Input Parameter:
2409: . v - the vector
2411: Output Parameter:
2412: . flg - flag indicating whether the boundtocpu flag will be propagated
2414: Level: developer
2416: .seealso: [](ch_vectors), `Vec`, `VecSetBindingPropagates()`
2417: @*/
2418: PetscErrorCode VecGetBindingPropagates(Vec v, PetscBool *flg)
2419: {
2420: PetscFunctionBegin;
2422: PetscAssertPointer(flg, 2);
2423: #if PetscDefined(HAVE_VIENNACL) || PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
2424: *flg = v->bindingpropagates;
2425: #else
2426: *flg = PETSC_FALSE;
2427: #endif
2428: PetscFunctionReturn(PETSC_SUCCESS);
2429: }
2431: /*@
2432: VecSetPinnedMemoryMin - Set the minimum data size for which pinned memory will be used for host (CPU) allocations.
2434: Logically Collective
2436: Input Parameters:
2437: + v - the vector
2438: - mbytes - minimum data size in bytes
2440: Options Database Key:
2441: . -vec_pinned_memory_min size - minimum size (in bytes) for an allocation to use pinned memory on host.
2443: Level: developer
2445: Note:
2446: Specifying -1 ensures that pinned memory will never be used.
2448: .seealso: [](ch_vectors), `Vec`, `VecGetPinnedMemoryMin()`
2449: @*/
2450: PetscErrorCode VecSetPinnedMemoryMin(Vec v, size_t mbytes)
2451: {
2452: PetscFunctionBegin;
2454: #if PetscDefined(HAVE_DEVICE)
2455: v->minimum_bytes_pinned_memory = mbytes;
2456: #endif
2457: PetscFunctionReturn(PETSC_SUCCESS);
2458: }
2460: /*@
2461: VecGetPinnedMemoryMin - Get the minimum data size for which pinned memory will be used for host (CPU) allocations.
2463: Logically Collective
2465: Input Parameter:
2466: . v - the vector
2468: Output Parameter:
2469: . mbytes - minimum data size in bytes
2471: Level: developer
2473: .seealso: [](ch_vectors), `Vec`, `VecSetPinnedMemoryMin()`
2474: @*/
2475: PetscErrorCode VecGetPinnedMemoryMin(Vec v, size_t *mbytes)
2476: {
2477: PetscFunctionBegin;
2479: PetscAssertPointer(mbytes, 2);
2480: #if PetscDefined(HAVE_DEVICE)
2481: *mbytes = v->minimum_bytes_pinned_memory;
2482: #endif
2483: PetscFunctionReturn(PETSC_SUCCESS);
2484: }
2486: /*@
2487: VecGetOffloadMask - Get the offload mask of a `Vec`
2489: Not Collective
2491: Input Parameter:
2492: . v - the vector
2494: Output Parameter:
2495: . mask - corresponding `PetscOffloadMask` enum value.
2497: Level: intermediate
2499: .seealso: [](ch_vectors), `Vec`, `VecCreateSeqCUDA()`, `VecCreateSeqViennaCL()`, `VecGetArray()`, `VecGetType()`
2500: @*/
2501: PetscErrorCode VecGetOffloadMask(Vec v, PetscOffloadMask *mask)
2502: {
2503: PetscFunctionBegin;
2505: PetscAssertPointer(mask, 2);
2506: *mask = v->offloadmask;
2507: PetscFunctionReturn(PETSC_SUCCESS);
2508: }
2510: #if !PetscDefined(HAVE_VIENNACL)
2511: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLContext(Vec v, PETSC_UINTPTR_T *ctx)
2512: {
2513: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_context");
2514: }
2516: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLQueue(Vec v, PETSC_UINTPTR_T *queue)
2517: {
2518: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_command_queue");
2519: }
2521: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMem(Vec v, PETSC_UINTPTR_T *queue)
2522: {
2523: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2524: }
2526: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMemRead(Vec v, PETSC_UINTPTR_T *queue)
2527: {
2528: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2529: }
2531: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMemWrite(Vec v, PETSC_UINTPTR_T *queue)
2532: {
2533: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2534: }
2536: PETSC_EXTERN PetscErrorCode VecViennaCLRestoreCLMemWrite(Vec v)
2537: {
2538: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to restore a Vec's cl_mem");
2539: }
2540: #endif
2542: 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)
2543: {
2544: const PetscScalar *u, *y;
2545: const PetscScalar *atola = NULL, *rtola = NULL, *erra = NULL;
2546: PetscInt n, n_loc = 0, na_loc = 0, nr_loc = 0;
2547: PetscReal nrm = 0, nrma = 0, nrmr = 0, err_loc[6];
2549: PetscFunctionBegin;
2550: #define SkipSmallValue(a, b, tol) \
2551: if (PetscAbsScalar(a) < tol || PetscAbsScalar(b) < tol) continue
2553: PetscCall(VecGetLocalSize(U, &n));
2554: PetscCall(VecGetArrayRead(U, &u));
2555: PetscCall(VecGetArrayRead(Y, &y));
2556: if (E) PetscCall(VecGetArrayRead(E, &erra));
2557: if (vatol) PetscCall(VecGetArrayRead(vatol, &atola));
2558: if (vrtol) PetscCall(VecGetArrayRead(vrtol, &rtola));
2559: for (PetscInt i = 0; i < n; i++) {
2560: PetscReal err, tol, tola, tolr;
2562: SkipSmallValue(y[i], u[i], ignore_max);
2563: atol = atola ? PetscRealPart(atola[i]) : atol;
2564: rtol = rtola ? PetscRealPart(rtola[i]) : rtol;
2565: err = erra ? PetscAbsScalar(erra[i]) : PetscAbsScalar(y[i] - u[i]);
2566: tola = atol;
2567: tolr = rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
2568: tol = tola + tolr;
2569: if (tola > 0.) {
2570: if (wnormtype == NORM_INFINITY) nrma = PetscMax(nrma, err / tola);
2571: else nrma += PetscSqr(err / tola);
2572: na_loc++;
2573: }
2574: if (tolr > 0.) {
2575: if (wnormtype == NORM_INFINITY) nrmr = PetscMax(nrmr, err / tolr);
2576: else nrmr += PetscSqr(err / tolr);
2577: nr_loc++;
2578: }
2579: if (tol > 0.) {
2580: if (wnormtype == NORM_INFINITY) nrm = PetscMax(nrm, err / tol);
2581: else nrm += PetscSqr(err / tol);
2582: n_loc++;
2583: }
2584: }
2585: if (E) PetscCall(VecRestoreArrayRead(E, &erra));
2586: if (vatol) PetscCall(VecRestoreArrayRead(vatol, &atola));
2587: if (vrtol) PetscCall(VecRestoreArrayRead(vrtol, &rtola));
2588: PetscCall(VecRestoreArrayRead(U, &u));
2589: PetscCall(VecRestoreArrayRead(Y, &y));
2590: #undef SkipSmallValue
2592: err_loc[0] = nrm;
2593: err_loc[1] = nrma;
2594: err_loc[2] = nrmr;
2595: err_loc[3] = (PetscReal)n_loc;
2596: err_loc[4] = (PetscReal)na_loc;
2597: err_loc[5] = (PetscReal)nr_loc;
2598: if (wnormtype == NORM_2) {
2599: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, err_loc, 6, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)U)));
2600: } else {
2601: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, err_loc, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)U)));
2602: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, err_loc + 3, 3, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)U)));
2603: }
2604: if (wnormtype == NORM_2) {
2605: *norm = PetscSqrtReal(err_loc[0]);
2606: *norma = PetscSqrtReal(err_loc[1]);
2607: *normr = PetscSqrtReal(err_loc[2]);
2608: } else {
2609: *norm = err_loc[0];
2610: *norma = err_loc[1];
2611: *normr = err_loc[2];
2612: }
2613: *norm_loc = (PetscInt)err_loc[3];
2614: *norma_loc = (PetscInt)err_loc[4];
2615: *normr_loc = (PetscInt)err_loc[5];
2616: PetscFunctionReturn(PETSC_SUCCESS);
2617: }
2619: /*@
2620: VecErrorWeightedNorms - compute a weighted norm of the difference between two vectors
2622: Collective
2624: Input Parameters:
2625: + U - first vector to be compared
2626: . Y - second vector to be compared
2627: . E - optional third vector representing the error (if not provided, the error is ||U-Y||)
2628: . wnormtype - norm type
2629: . atol - scalar for absolute tolerance
2630: . vatol - vector representing per-entry absolute tolerances (can be ``NULL``)
2631: . rtol - scalar for relative tolerance
2632: . vrtol - vector representing per-entry relative tolerances (can be ``NULL``)
2633: - ignore_max - ignore values smaller than this value in absolute terms.
2635: Output Parameters:
2636: + norm - weighted norm
2637: . norm_loc - number of vector locations used for the weighted norm
2638: . norma - weighted norm based on the absolute tolerance
2639: . norma_loc - number of vector locations used for the absolute weighted norm
2640: . normr - weighted norm based on the relative tolerance
2641: - normr_loc - number of vector locations used for the relative weighted norm
2643: Level: developer
2645: Notes:
2646: This is primarily used for computing weighted local truncation errors in ``TS``.
2648: .seealso: [](ch_vectors), `Vec`, `NormType`, `TSErrorWeightedNorm()`, `TSErrorWeightedENorm()`
2649: @*/
2650: 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)
2651: {
2652: PetscFunctionBegin;
2657: if (E) {
2660: }
2663: if (vatol) {
2666: }
2668: if (vrtol) {
2671: }
2673: PetscAssertPointer(norm, 10);
2674: PetscAssertPointer(norm_loc, 11);
2675: PetscAssertPointer(norma, 12);
2676: PetscAssertPointer(norma_loc, 13);
2677: PetscAssertPointer(normr, 14);
2678: PetscAssertPointer(normr_loc, 15);
2679: PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)U), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
2681: /* There are potentially 5 vectors involved, some of them may happen to be of different type or bound to cpu.
2682: Here we check that they all implement the same operation and call it if so.
2683: Otherwise, we call the _Basic implementation that always works (provided VecGetArrayRead is implemented). */
2684: PetscBool sameop = (PetscBool)(U->ops->errorwnorm && U->ops->errorwnorm == Y->ops->errorwnorm);
2685: if (sameop && E) sameop = (PetscBool)(U->ops->errorwnorm == E->ops->errorwnorm);
2686: if (sameop && vatol) sameop = (PetscBool)(U->ops->errorwnorm == vatol->ops->errorwnorm);
2687: if (sameop && vrtol) sameop = (PetscBool)(U->ops->errorwnorm == vrtol->ops->errorwnorm);
2688: if (sameop) PetscUseTypeMethod(U, errorwnorm, Y, E, wnormtype, atol, vatol, rtol, vrtol, ignore_max, norm, norm_loc, norma, norma_loc, normr, normr_loc);
2689: else PetscCall(VecErrorWeightedNorms_Basic(U, Y, E, wnormtype, atol, vatol, rtol, vrtol, ignore_max, norm, norm_loc, norma, norma_loc, normr, normr_loc));
2690: PetscFunctionReturn(PETSC_SUCCESS);
2691: }