Actual source code: vector.c
1: /*
2: Provides the interface functions for vector operations that do NOT have PetscScalar/PetscReal in the signature
3: These are the vector functions the user calls.
4: */
5: #include <petsc/private/vecimpl.h>
6: #include <petsc/private/deviceimpl.h>
8: /* Logging support */
9: PetscClassId VEC_CLASSID;
10: PetscLogEvent VEC_View, VEC_Max, VEC_Min, VEC_Dot, VEC_MDot, VEC_TDot;
11: PetscLogEvent VEC_Norm, VEC_Normalize, VEC_Scale, VEC_Shift, VEC_Copy, VEC_Set, VEC_AXPY, VEC_AYPX, VEC_WAXPY;
12: PetscLogEvent VEC_MTDot, VEC_MAXPY, VEC_Swap, VEC_AssemblyBegin, VEC_ScatterBegin, VEC_ScatterEnd;
13: PetscLogEvent VEC_AssemblyEnd, VEC_PointwiseMult, VEC_PointwiseDivide, VEC_Reciprocal, VEC_SetValues, VEC_Load, VEC_SetPreallocateCOO, VEC_SetValuesCOO;
14: PetscLogEvent VEC_SetRandom, VEC_ReduceArithmetic, VEC_ReduceCommunication, VEC_ReduceBegin, VEC_ReduceEnd, VEC_Ops;
15: PetscLogEvent VEC_DotNorm2, VEC_AXPBYPCZ;
16: PetscLogEvent VEC_ViennaCLCopyFromGPU, VEC_ViennaCLCopyToGPU;
17: PetscLogEvent VEC_CUDACopyFromGPU, VEC_CUDACopyToGPU;
18: PetscLogEvent VEC_HIPCopyFromGPU, VEC_HIPCopyToGPU;
20: /*@
21: VecStashGetInfo - Gets how many values are currently in the vector stash, i.e. need
22: to be communicated to other processors during the `VecAssemblyBegin()`/`VecAssemblyEnd()` process
24: Not Collective
26: Input Parameter:
27: . vec - the vector
29: Output Parameters:
30: + nstash - the size of the stash
31: . reallocs - the number of additional mallocs incurred in building the stash
32: . bnstash - the size of the block stash
33: - breallocs - the number of additional mallocs incurred in building the block stash (from `VecSetValuesBlocked()`)
35: Level: advanced
37: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecStashSetInitialSize()`, `VecStashView()`
38: @*/
39: PetscErrorCode VecStashGetInfo(Vec vec, PetscInt *nstash, PetscInt *reallocs, PetscInt *bnstash, PetscInt *breallocs)
40: {
41: PetscFunctionBegin;
42: PetscCall(VecStashGetInfo_Private(&vec->stash, nstash, reallocs));
43: PetscCall(VecStashGetInfo_Private(&vec->bstash, bnstash, breallocs));
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: /*@
48: VecSetLocalToGlobalMapping - Sets a local numbering to global numbering used
49: by the routine `VecSetValuesLocal()` to allow users to insert vector entries
50: using a local (per-processor) numbering.
52: Logically Collective
54: Input Parameters:
55: + x - vector
56: - mapping - mapping created with `ISLocalToGlobalMappingCreate()` or `ISLocalToGlobalMappingCreateIS()`
58: Level: intermediate
60: Notes:
61: All vectors obtained with `VecDuplicate()` from this vector inherit the same mapping.
63: Vectors obtained with `DMCreateGlobaVector()` will often have this attribute attached to the vector so this call is not needed
65: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesLocal()`,
66: `VecGetLocalToGlobalMapping()`, `VecSetValuesBlockedLocal()`
67: @*/
68: PetscErrorCode VecSetLocalToGlobalMapping(Vec x, ISLocalToGlobalMapping mapping)
69: {
70: PetscFunctionBegin;
73: if (x->ops->setlocaltoglobalmapping) PetscUseTypeMethod(x, setlocaltoglobalmapping, mapping);
74: else PetscCall(PetscLayoutSetISLocalToGlobalMapping(x->map, mapping));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: /*@
79: VecGetLocalToGlobalMapping - Gets the local-to-global numbering set by `VecSetLocalToGlobalMapping()`
81: Not Collective
83: Input Parameter:
84: . X - the vector
86: Output Parameter:
87: . mapping - the mapping
89: Level: advanced
91: .seealso: [](ch_vectors), `Vec`, `VecSetValuesLocal()`, `VecSetLocalToGlobalMapping()`
92: @*/
93: PetscErrorCode VecGetLocalToGlobalMapping(Vec X, ISLocalToGlobalMapping *mapping)
94: {
95: PetscFunctionBegin;
98: PetscAssertPointer(mapping, 2);
99: if (X->ops->getlocaltoglobalmapping) PetscUseTypeMethod(X, getlocaltoglobalmapping, mapping);
100: else *mapping = X->map->mapping;
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: /*@
105: VecAssemblyBegin - Begins assembling the vector; that is ensuring all the vector's entries are stored on the correct MPI process. This routine should
106: be called after completing all calls to `VecSetValues()`.
108: Collective
110: Input Parameter:
111: . vec - the vector
113: Level: beginner
115: .seealso: [](ch_vectors), `Vec`, `VecAssemblyEnd()`, `VecSetValues()`
116: @*/
117: PetscErrorCode VecAssemblyBegin(Vec vec)
118: {
119: PetscFunctionBegin;
122: PetscCall(VecStashViewFromOptions(vec, NULL, "-vec_view_stash"));
123: PetscCall(PetscLogEventBegin(VEC_AssemblyBegin, vec, 0, 0, 0));
124: PetscTryTypeMethod(vec, assemblybegin);
125: PetscCall(PetscLogEventEnd(VEC_AssemblyBegin, vec, 0, 0, 0));
126: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: /*@
131: VecAssemblyEnd - Completes assembling the vector. This routine should be called after `VecAssemblyBegin()`.
133: Collective
135: Input Parameter:
136: . vec - the vector
138: Options Database Keys:
139: + -vec_view - Prints vector in `PETSC_VIEWER_DEFAULT` format
140: . -vec_view ::ascii_matlab - Prints vector in `PETSC_VIEWER_ASCII_MATLAB` format to stdout
141: . -vec_view matlab:filename - Prints vector in MATLAB .mat file to filename (requires PETSc configured with --with-matlab)
142: . -vec_view draw - Activates vector viewing using drawing tools
143: . -display name - Sets display name (default is host)
144: . -draw_pause sec - Sets number of seconds to pause after display
145: - -vec_view socket - Activates vector viewing using a socket
147: Level: beginner
149: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecSetValues()`
150: @*/
151: PetscErrorCode VecAssemblyEnd(Vec vec)
152: {
153: PetscFunctionBegin;
155: PetscCall(PetscLogEventBegin(VEC_AssemblyEnd, vec, 0, 0, 0));
157: PetscTryTypeMethod(vec, assemblyend);
158: PetscCall(PetscLogEventEnd(VEC_AssemblyEnd, vec, 0, 0, 0));
159: PetscCall(VecViewFromOptions(vec, NULL, "-vec_view"));
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: /*@
164: VecSetPreallocationCOO - set preallocation for a vector using a coordinate format of the entries with global indices
166: Collective
168: Input Parameters:
169: + x - vector being preallocated
170: . ncoo - number of entries
171: - coo_i - entry indices
173: Level: beginner
175: Notes:
176: This and `VecSetValuesCOO()` provide an alternative API to using `VecSetValues()` to provide vector values.
178: This API is particularly efficient for use on GPUs.
180: Entries can be repeated, see `VecSetValuesCOO()`. Negative indices are not allowed unless vector option `VEC_IGNORE_NEGATIVE_INDICES` is set,
181: in which case they, along with the corresponding entries in `VecSetValuesCOO()`, are ignored. If vector option `VEC_NO_OFF_PROC_ENTRIES` is set,
182: remote entries are ignored, otherwise, they will be properly added or inserted to the vector.
184: The array coo_i[] may be freed immediately after calling this function.
186: .seealso: [](ch_vectors), `Vec`, `VecSetValuesCOO()`, `VecSetPreallocationCOOLocal()`
187: @*/
188: PetscErrorCode VecSetPreallocationCOO(Vec x, PetscCount ncoo, const PetscInt coo_i[])
189: {
190: PetscFunctionBegin;
193: if (ncoo) PetscAssertPointer(coo_i, 3);
194: PetscCall(PetscLogEventBegin(VEC_SetPreallocateCOO, x, 0, 0, 0));
195: PetscCall(PetscLayoutSetUp(x->map));
196: if (x->ops->setpreallocationcoo) {
197: PetscUseTypeMethod(x, setpreallocationcoo, ncoo, coo_i);
198: } else {
199: PetscInt ncoo_i;
200: IS is_coo_i;
202: PetscCall(PetscIntCast(ncoo, &ncoo_i));
203: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncoo_i, coo_i, PETSC_COPY_VALUES, &is_coo_i));
204: PetscCall(PetscObjectCompose((PetscObject)x, "__PETSc_coo_i", (PetscObject)is_coo_i));
205: PetscCall(ISDestroy(&is_coo_i));
206: }
207: PetscCall(PetscLogEventEnd(VEC_SetPreallocateCOO, x, 0, 0, 0));
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: /*@
212: VecSetPreallocationCOOLocal - set preallocation for vectors using a coordinate format of the entries with local indices
214: Collective
216: Input Parameters:
217: + x - vector being preallocated
218: . ncoo - number of entries
219: - coo_i - row indices (local numbering; may be modified)
221: Level: beginner
223: Notes:
224: This and `VecSetValuesCOO()` provide an alternative API to using `VecSetValuesLocal()` to provide vector values.
226: This API is particularly efficient for use on GPUs.
228: The local indices are translated using the local to global mapping, thus `VecSetLocalToGlobalMapping()` must have been
229: called prior to this function.
231: The indices coo_i may be modified within this function. They might be translated to corresponding global
232: indices, but the caller should not rely on them having any specific value after this function returns. The arrays
233: can be freed or reused immediately after this function returns.
235: Entries can be repeated. Negative indices and remote indices might be allowed. see `VecSetPreallocationCOO()`.
237: .seealso: [](ch_vectors), `Vec`, `VecSetPreallocationCOO()`, `VecSetValuesCOO()`
238: @*/
239: PetscErrorCode VecSetPreallocationCOOLocal(Vec x, PetscCount ncoo, PetscInt coo_i[])
240: {
241: PetscInt ncoo_i;
242: ISLocalToGlobalMapping ltog;
244: PetscFunctionBegin;
247: if (ncoo) PetscAssertPointer(coo_i, 3);
248: PetscCall(PetscIntCast(ncoo, &ncoo_i));
249: PetscCall(PetscLayoutSetUp(x->map));
250: PetscCall(VecGetLocalToGlobalMapping(x, <og));
251: if (ltog) PetscCall(ISLocalToGlobalMappingApply(ltog, ncoo_i, coo_i, coo_i));
252: PetscCall(VecSetPreallocationCOO(x, ncoo, coo_i));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: /*@
257: VecSetValuesCOO - set values at once in a vector preallocated using `VecSetPreallocationCOO()`
259: Collective
261: Input Parameters:
262: + x - vector being set
263: . coo_v - the value array
264: - imode - the insert mode
266: Level: beginner
268: Note:
269: This and `VecSetPreallocationCOO() or ``VecSetPreallocationCOOLocal()` provide an alternative API to using `VecSetValues()` to provide vector values.
271: This API is particularly efficient for use on GPUs.
273: The values must follow the order of the indices prescribed with `VecSetPreallocationCOO()` or `VecSetPreallocationCOOLocal()`.
274: When repeated entries are specified in the COO indices the `coo_v` values are first properly summed, regardless of the value of `imode`.
275: The imode flag indicates if `coo_v` must be added to the current values of the vector (`ADD_VALUES`) or overwritten (`INSERT_VALUES`).
276: `VecAssemblyBegin()` and `VecAssemblyEnd()` do not need to be called after this routine. It automatically handles the assembly process.
278: .seealso: [](ch_vectors), `Vec`, `VecSetPreallocationCOO()`, `VecSetPreallocationCOOLocal()`, `VecSetValues()`
279: @*/
280: PetscErrorCode VecSetValuesCOO(Vec x, const PetscScalar coo_v[], InsertMode imode)
281: {
282: PetscFunctionBegin;
286: PetscCall(PetscLogEventBegin(VEC_SetValuesCOO, x, 0, 0, 0));
287: if (x->ops->setvaluescoo) {
288: PetscUseTypeMethod(x, setvaluescoo, coo_v, imode);
289: PetscCall(PetscObjectStateIncrease((PetscObject)x));
290: } else {
291: IS is_coo_i;
292: const PetscInt *coo_i;
293: PetscInt ncoo;
294: PetscMemType mtype;
296: PetscCall(PetscGetMemType(coo_v, &mtype));
297: PetscCheck(mtype == PETSC_MEMTYPE_HOST, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_WRONG, "The basic VecSetValuesCOO() only supports v[] on host");
298: PetscCall(PetscObjectQuery((PetscObject)x, "__PETSc_coo_i", (PetscObject *)&is_coo_i));
299: PetscCheck(is_coo_i, PetscObjectComm((PetscObject)x), PETSC_ERR_COR, "Missing coo_i IS");
300: PetscCall(ISGetLocalSize(is_coo_i, &ncoo));
301: PetscCall(ISGetIndices(is_coo_i, &coo_i));
302: if (imode != ADD_VALUES) PetscCall(VecZeroEntries(x));
303: PetscCall(VecSetValues(x, ncoo, coo_i, coo_v, ADD_VALUES));
304: PetscCall(ISRestoreIndices(is_coo_i, &coo_i));
305: PetscCall(VecAssemblyBegin(x));
306: PetscCall(VecAssemblyEnd(x));
307: }
308: PetscCall(PetscLogEventEnd(VEC_SetValuesCOO, x, 0, 0, 0));
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: 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))
313: {
314: PetscErrorCode (*async_fn)(Vec, Vec, Vec, PetscDeviceContext) = NULL;
316: PetscFunctionBegin;
323: PetscCheckSameTypeAndComm(x, 2, y, 3);
324: PetscCheckSameTypeAndComm(y, 3, w, 1);
325: VecCheckSameSize(w, 1, x, 2);
326: VecCheckSameSize(w, 1, y, 3);
327: VecCheckAssembled(x);
328: VecCheckAssembled(y);
329: PetscCall(VecSetErrorIfLocked(w, 1));
332: if (dctx) PetscCall(PetscObjectQueryFunction((PetscObject)w, async_name, &async_fn));
333: if (event) PetscCall(PetscLogEventBegin(event, x, y, w, 0));
334: if (async_fn) PetscCall((*async_fn)(w, x, y, dctx));
335: else PetscCall((*pointwise_op)(w, x, y));
336: if (event) PetscCall(PetscLogEventEnd(event, x, y, w, 0));
337: PetscCall(PetscObjectStateIncrease((PetscObject)w));
338: PetscFunctionReturn(PETSC_SUCCESS);
339: }
341: PetscErrorCode VecPointwiseMaxAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
342: {
343: PetscFunctionBegin;
344: // REVIEW ME: no log event?
345: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMax), w->ops->pointwisemax));
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: /*@
350: VecPointwiseMax - Computes the component-wise maximum `w[i] = max(x[i], y[i])`.
352: Logically Collective
354: Input Parameters:
355: + x - the first input vector
356: - y - the second input vector
358: Output Parameter:
359: . w - the result
361: Level: advanced
363: Notes:
364: Any subset of the `x`, `y`, and `w` may be the same vector.
366: For complex numbers compares only the real part
368: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
369: @*/
370: PetscErrorCode VecPointwiseMax(Vec w, Vec x, Vec y)
371: {
372: PetscFunctionBegin;
373: PetscCall(VecPointwiseMaxAsync_Private(w, x, y, NULL));
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }
377: PetscErrorCode VecPointwiseMinAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
378: {
379: PetscFunctionBegin;
380: // REVIEW ME: no log event?
381: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMin), w->ops->pointwisemin));
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: /*@
386: VecPointwiseMin - Computes the component-wise minimum `w[i] = min(x[i], y[i])`.
388: Logically Collective
390: Input Parameters:
391: + x - the first input vector
392: - y - the second input vector
394: Output Parameter:
395: . w - the result
397: Level: advanced
399: Notes:
400: Any subset of the `x`, `y`, and `w` may be the same vector.
402: For complex numbers compares only the real part
404: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
405: @*/
406: PetscErrorCode VecPointwiseMin(Vec w, Vec x, Vec y)
407: {
408: PetscFunctionBegin;
409: PetscCall(VecPointwiseMinAsync_Private(w, x, y, NULL));
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: PetscErrorCode VecPointwiseMaxAbsAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
414: {
415: PetscFunctionBegin;
416: // REVIEW ME: no log event?
417: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMaxAbs), w->ops->pointwisemaxabs));
418: PetscFunctionReturn(PETSC_SUCCESS);
419: }
421: /*@
422: VecPointwiseMaxAbs - Computes the component-wise maximum of the absolute values `w[i] = max(abs(x[i]), abs(y[i]))`.
424: Logically Collective
426: Input Parameters:
427: + x - the first input vector
428: - y - the second input vector
430: Output Parameter:
431: . w - the result
433: Level: advanced
435: Notes:
436: Any subset of the `x`, `y`, and `w` may be the same vector.
438: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMin()`, `VecPointwiseMax()`, `VecMaxPointwiseDivide()`
439: @*/
440: PetscErrorCode VecPointwiseMaxAbs(Vec w, Vec x, Vec y)
441: {
442: PetscFunctionBegin;
443: PetscCall(VecPointwiseMaxAbsAsync_Private(w, x, y, NULL));
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: PetscErrorCode VecPointwiseDivideAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
448: {
449: PetscFunctionBegin;
450: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, VEC_PointwiseDivide, VecAsyncFnName(PointwiseDivide), w->ops->pointwisedivide));
451: PetscFunctionReturn(PETSC_SUCCESS);
452: }
454: /*@
455: VecPointwiseDivide - Computes the component-wise division `w[i] = x[i] / y[i]`.
457: Logically Collective
459: Input Parameters:
460: + x - the numerator vector
461: - y - the denominator vector
463: Output Parameter:
464: . w - the result
466: Level: advanced
468: Note:
469: Any subset of the `x`, `y`, and `w` may be the same vector.
471: .seealso: [](ch_vectors), `Vec`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
472: @*/
473: PetscErrorCode VecPointwiseDivide(Vec w, Vec x, Vec y)
474: {
475: PetscFunctionBegin;
476: PetscCall(VecPointwiseDivideAsync_Private(w, x, y, NULL));
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }
480: #define VEC_POINTWISE_SIGN_LOOP(y, x, n, func) \
481: PetscPragmaSIMD \
482: for (PetscInt i = 0; i < (n); i++) (y)[i] = func(PetscRealPart((x)[i]))
484: #define VEC_POINTWISE_SIGN_DISPATCH(y, x, n, sign_type) \
485: do { \
486: switch (sign_type) { \
487: case VEC_SIGN_ZERO_TO_ZERO: \
488: VEC_POINTWISE_SIGN_LOOP(y, x, n, VecSignZeroToZero_Private); \
489: break; \
490: case VEC_SIGN_ZERO_TO_SIGNED_ZERO: \
491: VEC_POINTWISE_SIGN_LOOP(y, x, n, VecSignZeroToSignedZero_Private); \
492: break; \
493: case VEC_SIGN_ZERO_TO_SIGNED_UNIT: \
494: VEC_POINTWISE_SIGN_LOOP(y, x, n, VecSignZeroToSignedUnit_Private); \
495: break; \
496: default: \
497: PetscUnreachable(); \
498: } \
499: } while (0)
501: PetscErrorCode VecPointwiseSignAsync_Private(Vec y, Vec x, VecSignMode sign_type, PetscDeviceContext dctx)
502: {
503: PetscOffloadMask mask;
504: PetscBool is_host;
505: PetscErrorCode (*async_fn)(Vec, Vec, VecSignMode, PetscDeviceContext) = NULL;
507: PetscFunctionBegin;
512: VecCheckSameSize(y, 1, x, 2);
513: VecCheckAssembled(x);
514: VecCheckAssembled(y);
515: PetscCall(VecSetErrorIfLocked(y, 1));
517: PetscCall(VecGetOffloadMask(x, &mask));
518: is_host = PetscOffloadHost(mask) ? PETSC_TRUE : PETSC_FALSE;
519: if (!is_host) PetscCall(PetscObjectQueryFunction((PetscObject)y, VEC_ASYNC_FN_NAME("PointwiseSign"), &async_fn));
520: if (async_fn) PetscCall((*async_fn)(y, x, sign_type, dctx));
521: else {
522: PetscInt n;
524: PetscCall(VecGetLocalSize(y, &n));
525: if (y == x) {
526: PetscScalar *_y;
528: PetscCall(VecGetArray(y, &_y));
529: VEC_POINTWISE_SIGN_DISPATCH(_y, _y, n, sign_type);
530: PetscCall(VecRestoreArray(y, &_y));
531: } else {
532: PetscScalar *_y;
533: const PetscScalar *_x;
535: PetscCall(VecGetArrayWrite(y, &_y));
536: PetscCall(VecGetArrayRead(x, &_x));
537: VEC_POINTWISE_SIGN_DISPATCH(_y, _x, n, sign_type);
538: PetscCall(VecRestoreArrayRead(x, &_x));
539: PetscCall(VecRestoreArrayWrite(y, &_y));
540: }
541: }
542: PetscCall(PetscObjectStateIncrease((PetscObject)y));
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: /*@
547: VecPointwiseSign - Computes the component-wise sign `y[i] = sign(x[i])`.
549: Logically Collective
551: Input Parameters:
552: + x - the input vector
553: - sign_type - `VecSignMode` indicating how the function should map zero values.
555: Output Parameter:
556: . y - the sign vector of `x`
558: Level: beginner
560: .seealso: [](ch_vectors), `Vec`, `VecSignMode`
561: @*/
562: PetscErrorCode VecPointwiseSign(Vec y, Vec x, VecSignMode sign_type)
563: {
564: PetscFunctionBegin;
565: PetscCall(VecPointwiseSignAsync_Private(y, x, sign_type, NULL));
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: PetscErrorCode VecPointwiseMultAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
570: {
571: PetscFunctionBegin;
573: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, VEC_PointwiseMult, VecAsyncFnName(PointwiseMult), w->ops->pointwisemult));
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: /*@
578: VecPointwiseMult - Computes the component-wise multiplication `w[i] = x[i] * y[i]`.
580: Logically Collective
582: Input Parameters:
583: + x - the first vector
584: - y - the second vector
586: Output Parameter:
587: . w - the result
589: Level: advanced
591: Note:
592: Any subset of the `x`, `y`, and `w` may be the same vector.
594: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
595: @*/
596: PetscErrorCode VecPointwiseMult(Vec w, Vec x, Vec y)
597: {
598: PetscFunctionBegin;
599: PetscCall(VecPointwiseMultAsync_Private(w, x, y, NULL));
600: PetscFunctionReturn(PETSC_SUCCESS);
601: }
603: /*@
604: VecDuplicate - Creates a new vector of the same type as an existing vector.
606: Collective
608: Input Parameter:
609: . v - a vector to mimic
611: Output Parameter:
612: . newv - location to put new vector
614: Level: beginner
616: Notes:
617: `VecDuplicate()` DOES NOT COPY the vector entries, but rather allocates storage
618: for the new vector. Use `VecCopy()` to copy a vector.
620: Use `VecDestroy()` to free the space. Use `VecDuplicateVecs()` to get several
621: vectors.
623: .seealso: [](ch_vectors), `Vec`, `VecDestroy()`, `VecDuplicateVecs()`, `VecCreate()`, `VecCopy()`
624: @*/
625: PetscErrorCode VecDuplicate(Vec v, Vec *newv)
626: {
627: PetscFunctionBegin;
629: PetscAssertPointer(newv, 2);
631: PetscUseTypeMethod(v, duplicate, newv);
632: #if PetscDefined(HAVE_DEVICE)
633: if (v->boundtocpu && v->bindingpropagates) {
634: PetscCall(VecSetBindingPropagates(*newv, PETSC_TRUE));
635: PetscCall(VecBindToCPU(*newv, PETSC_TRUE));
636: }
637: #endif
638: PetscCall(PetscObjectStateIncrease((PetscObject)*newv));
639: PetscFunctionReturn(PETSC_SUCCESS);
640: }
642: /*@
643: VecDestroy - Destroys a vector.
645: Collective
647: Input Parameter:
648: . v - the vector
650: Level: beginner
652: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecDuplicate()`, `VecDestroyVecs()`
653: @*/
654: PetscErrorCode VecDestroy(Vec *v)
655: {
656: PetscFunctionBegin;
657: PetscAssertPointer(v, 1);
658: if (!*v) PetscFunctionReturn(PETSC_SUCCESS);
660: if (--((PetscObject)*v)->refct > 0) {
661: *v = NULL;
662: PetscFunctionReturn(PETSC_SUCCESS);
663: }
665: PetscCall(PetscObjectSAWsViewOff((PetscObject)*v));
666: /* destroy the internal part */
667: PetscTryTypeMethod(*v, destroy);
668: PetscCall(PetscFree((*v)->defaultrandtype));
669: /* destroy the external/common part */
670: PetscCall(PetscLayoutDestroy(&(*v)->map));
671: PetscCall(PetscHeaderDestroy(v));
672: PetscFunctionReturn(PETSC_SUCCESS);
673: }
675: /*@C
676: VecDuplicateVecs - Creates several vectors of the same type as an existing vector.
678: Collective
680: Input Parameters:
681: + m - the number of vectors to obtain
682: - v - a vector to mimic
684: Output Parameter:
685: . V - location to put pointer to array of vectors
687: Level: intermediate
689: Notes:
690: Use `VecDestroyVecs()` to free the space. Use `VecDuplicate()` to form a single
691: vector.
693: Some implementations ensure that the arrays accessed by each vector are contiguous in memory. Certain `VecMDot()` and `VecMAXPY()`
694: implementations utilize this property to use BLAS 2 operations for higher efficiency. This is especially useful in `KSPGMRES`, see
695: `KSPGMRESSetPreAllocateVectors()`.
697: Fortran Note:
698: .vb
699: Vec, pointer :: V(:)
700: .ve
702: .seealso: [](ch_vectors), `Vec`, [](ch_fortran), `VecDestroyVecs()`, `VecDuplicate()`, `VecCreate()`, `VecMDot()`, `VecMAXPY()`, `KSPGMRES`,
703: `KSPGMRESSetPreAllocateVectors()`
704: @*/
705: PetscErrorCode VecDuplicateVecs(Vec v, PetscInt m, Vec *V[])
706: {
707: PetscFunctionBegin;
709: PetscAssertPointer(V, 3);
711: PetscUseTypeMethod(v, duplicatevecs, m, V);
712: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
713: if (v->boundtocpu && v->bindingpropagates) {
714: PetscInt i;
716: for (i = 0; i < m; i++) {
717: /* Since ops->duplicatevecs might itself propagate the value of boundtocpu,
718: * avoid unnecessary overhead by only calling VecBindToCPU() if the vector isn't already bound. */
719: if (!(*V)[i]->boundtocpu) {
720: PetscCall(VecSetBindingPropagates((*V)[i], PETSC_TRUE));
721: PetscCall(VecBindToCPU((*V)[i], PETSC_TRUE));
722: }
723: }
724: }
725: #endif
726: PetscFunctionReturn(PETSC_SUCCESS);
727: }
729: /*@C
730: VecDestroyVecs - Frees a block of vectors obtained with `VecDuplicateVecs()`.
732: Collective
734: Input Parameters:
735: + m - the number of vectors previously obtained, if zero no vectors are destroyed
736: - vv - pointer to pointer to array of vector pointers, if `NULL` no vectors are destroyed
738: Level: intermediate
740: .seealso: [](ch_vectors), `Vec`, [](ch_fortran), `VecDuplicateVecs()`, `VecDestroyVecsf90()`
741: @*/
742: PetscErrorCode VecDestroyVecs(PetscInt m, Vec *vv[])
743: {
744: PetscFunctionBegin;
745: PetscAssertPointer(vv, 2);
746: PetscCheck(m >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Trying to destroy negative number of vectors %" PetscInt_FMT, m);
747: if (!m || !*vv) {
748: *vv = NULL;
749: PetscFunctionReturn(PETSC_SUCCESS);
750: }
753: PetscCall((*(**vv)->ops->destroyvecs)(m, *vv));
754: *vv = NULL;
755: PetscFunctionReturn(PETSC_SUCCESS);
756: }
758: /*@
759: VecViewFromOptions - View a vector based on values in the options database
761: Collective
763: Input Parameters:
764: + A - the vector
765: . obj - optional object that provides the options prefix for this viewing, use 'NULL' to use the prefix of `A`
766: - name - command line option
768: Level: intermediate
770: Note:
771: See `PetscObjectViewFromOptions()` to see the `PetscViewer` and PetscViewerFormat` available
773: .seealso: [](ch_vectors), `Vec`, `VecView`, `PetscObjectViewFromOptions()`, `VecCreate()`
774: @*/
775: PetscErrorCode VecViewFromOptions(Vec A, PeOp PetscObject obj, const char name[])
776: {
777: PetscFunctionBegin;
779: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
780: PetscFunctionReturn(PETSC_SUCCESS);
781: }
783: /*@
784: VecView - Views a vector object.
786: Collective
788: Input Parameters:
789: + vec - the vector
790: - viewer - an optional `PetscViewer` visualization context
792: Level: beginner
794: Notes:
795: The available visualization contexts include
796: + `PETSC_VIEWER_STDOUT_SELF` - for sequential vectors
797: . `PETSC_VIEWER_STDOUT_WORLD` - for parallel vectors created on `PETSC_COMM_WORLD`
798: - `PETSC_VIEWER_STDOUT`_(comm) - for parallel vectors created on MPI communicator comm
800: You can change the format the vector is printed using the
801: option `PetscViewerPushFormat()`.
803: The user can open alternative viewers with
804: + `PetscViewerASCIIOpen()` - Outputs vector to a specified file
805: . `PetscViewerBinaryOpen()` - Outputs vector in binary to a
806: specified file; corresponding input uses `VecLoad()`
807: . `PetscViewerDrawOpen()` - Outputs vector to an X window display
808: . `PetscViewerSocketOpen()` - Outputs vector to Socket viewer
809: - `PetscViewerHDF5Open()` - Outputs vector to HDF5 file viewer
811: The user can call `PetscViewerPushFormat()` to specify the output
812: format of ASCII printed objects (when using `PETSC_VIEWER_STDOUT_SELF`,
813: `PETSC_VIEWER_STDOUT_WORLD` and `PetscViewerASCIIOpen()`). Available formats include
814: + `PETSC_VIEWER_DEFAULT` - default, prints vector contents
815: . `PETSC_VIEWER_ASCII_MATLAB` - prints vector contents in MATLAB format
816: . `PETSC_VIEWER_ASCII_INDEX` - prints vector contents, including indices of vector elements
817: - `PETSC_VIEWER_ASCII_COMMON` - prints vector contents, using a
818: format common among all vector types
820: You can pass any number of vector objects, or other PETSc objects to the same viewer.
822: In the debugger you can do call `VecView`(v,0) to display the vector. (The same holds for any PETSc object viewer).
824: Notes for binary viewer:
825: If you pass multiple vectors to a binary viewer you can read them back in the same order
826: with `VecLoad()`.
828: If the blocksize of the vector is greater than one then you must provide a unique prefix to
829: the vector with `PetscObjectSetOptionsPrefix`((`PetscObject`)vec,"uniqueprefix"); BEFORE calling `VecView()` on the
830: vector to be stored and then set that same unique prefix on the vector that you pass to `VecLoad()`. The blocksize
831: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
832: filename. If you copy the binary file, make sure you copy the associated .info file with it.
834: See the manual page for `VecLoad()` on the exact format the binary viewer stores
835: the values in the file.
837: Notes for HDF5 Viewer:
838: The name of the `Vec` (given with `PetscObjectSetName()` is the name that is used
839: for the object in the HDF5 file. If you wish to store the same Vec into multiple
840: datasets in the same file (typically with different values), you must change its
841: name each time before calling the `VecView()`. To load the same vector,
842: the name of the Vec object passed to `VecLoad()` must be the same.
844: If the block size of the vector is greater than 1 then it is used as the first dimension in the HDF5 array.
845: If the function `PetscViewerHDF5SetBaseDimension2()`is called then even if the block size is one it will
846: be used as the first dimension in the HDF5 array (that is the HDF5 array will always be two dimensional)
847: See also `PetscViewerHDF5SetTimestep()` which adds an additional complication to reading and writing `Vec`
848: with the HDF5 viewer.
850: .seealso: [](ch_vectors), `Vec`, `VecViewFromOptions()`, `PetscViewerASCIIOpen()`, `PetscViewerDrawOpen()`, `PetscDrawLGCreate()`,
851: `PetscViewerSocketOpen()`, `PetscViewerBinaryOpen()`, `VecLoad()`, `PetscViewerCreate()`,
852: `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`, `PetscViewerHDF5SetTimestep()`
853: @*/
854: PetscErrorCode VecView(Vec vec, PetscViewer viewer)
855: {
856: PetscBool isascii;
857: PetscViewerFormat format;
858: PetscMPIInt size;
860: PetscFunctionBegin;
863: VecCheckAssembled(vec);
864: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec), &viewer));
866: PetscCall(PetscViewerGetFormat(viewer, &format));
867: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
868: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) PetscFunctionReturn(PETSC_SUCCESS);
870: PetscCheck(!vec->stash.n && !vec->bstash.n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call VecAssemblyBegin/End() before viewing this vector");
872: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
873: if (isascii) {
874: PetscInt rows, bs;
876: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)vec, viewer));
877: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
878: PetscCall(PetscViewerASCIIPushTab(viewer));
879: PetscCall(VecGetSize(vec, &rows));
880: PetscCall(VecGetBlockSize(vec, &bs));
881: if (bs != 1) {
882: PetscCall(PetscViewerASCIIPrintf(viewer, "length=%" PetscInt_FMT ", bs=%" PetscInt_FMT "\n", rows, bs));
883: } else {
884: PetscCall(PetscViewerASCIIPrintf(viewer, "length=%" PetscInt_FMT "\n", rows));
885: }
886: PetscCall(PetscViewerASCIIPopTab(viewer));
887: }
888: }
889: PetscCall(VecLockReadPush(vec));
890: PetscCall(PetscLogEventBegin(VEC_View, vec, viewer, 0, 0));
891: if ((format == PETSC_VIEWER_NATIVE || format == PETSC_VIEWER_LOAD_BALANCE) && vec->ops->viewnative) {
892: PetscUseTypeMethod(vec, viewnative, viewer);
893: } else {
894: PetscUseTypeMethod(vec, view, viewer);
895: }
896: PetscCall(VecLockReadPop(vec));
897: PetscCall(PetscLogEventEnd(VEC_View, vec, viewer, 0, 0));
898: PetscFunctionReturn(PETSC_SUCCESS);
899: }
901: #if defined(PETSC_USE_DEBUG)
902: #include <../src/sys/totalview/tv_data_display.h>
903: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
904: {
905: const PetscScalar *values;
906: char type[32];
908: TV_add_row("Local rows", "int", &v->map->n);
909: TV_add_row("Global rows", "int", &v->map->N);
910: TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
911: PetscCall(VecGetArrayRead((Vec)v, &values));
912: PetscCall(PetscSNPrintf(type, 32, "double[%" PetscInt_FMT "]", v->map->n));
913: TV_add_row("values", type, values);
914: PetscCall(VecRestoreArrayRead((Vec)v, &values));
915: return TV_format_OK;
916: }
917: #endif
919: /*@C
920: VecViewNative - Views a vector object with the original type specific viewer
922: Collective
924: Input Parameters:
925: + vec - the vector
926: - viewer - an optional `PetscViewer` visualization context
928: Level: developer
930: Note:
931: This can be used with, for example, vectors obtained with `DMCreateGlobalVector()` for a `DMDA` to display the vector
932: in the PETSc storage format (each MPI process values follow the previous MPI processes) instead of the "natural" grid
933: ordering.
935: .seealso: [](ch_vectors), `Vec`, `PetscViewerASCIIOpen()`, `PetscViewerDrawOpen()`, `PetscDrawLGCreate()`, `VecView()`,
936: `PetscViewerSocketOpen()`, `PetscViewerBinaryOpen()`, `VecLoad()`, `PetscViewerCreate()`,
937: `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`, `PetscViewerHDF5SetTimestep()`
938: @*/
939: PetscErrorCode VecViewNative(Vec vec, PetscViewer viewer)
940: {
941: PetscFunctionBegin;
944: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec), &viewer));
946: PetscUseTypeMethod(vec, viewnative, viewer);
947: PetscFunctionReturn(PETSC_SUCCESS);
948: }
950: /*@
951: VecGetSize - Returns the global number of elements of the vector.
953: Not Collective
955: Input Parameter:
956: . x - the vector
958: Output Parameter:
959: . size - the global length of the vector
961: Level: beginner
963: .seealso: [](ch_vectors), `Vec`, `VecGetLocalSize()`
964: @*/
965: PetscErrorCode VecGetSize(Vec x, PetscInt *size)
966: {
967: PetscFunctionBegin;
969: PetscAssertPointer(size, 2);
971: PetscUseTypeMethod(x, getsize, size);
972: PetscFunctionReturn(PETSC_SUCCESS);
973: }
975: /*@
976: VecGetLocalSize - Returns the number of elements of the vector stored
977: in local memory (that is on this MPI process)
979: Not Collective
981: Input Parameter:
982: . x - the vector
984: Output Parameter:
985: . size - the length of the local piece of the vector
987: Level: beginner
989: .seealso: [](ch_vectors), `Vec`, `VecGetSize()`
990: @*/
991: PetscErrorCode VecGetLocalSize(Vec x, PetscInt *size)
992: {
993: PetscFunctionBegin;
995: PetscAssertPointer(size, 2);
997: PetscUseTypeMethod(x, getlocalsize, size);
998: PetscFunctionReturn(PETSC_SUCCESS);
999: }
1001: /*@
1002: VecGetOwnershipRange - Returns the range of indices owned by
1003: this process. The vector is laid out with the
1004: first `n1` elements on the first processor, next `n2` elements on the
1005: second, etc. For certain parallel layouts this range may not be
1006: well defined.
1008: Not Collective
1010: Input Parameter:
1011: . x - the vector
1013: Output Parameters:
1014: + low - the first local element, pass in `NULL` if not interested
1015: - high - one more than the last local element, pass in `NULL` if not interested
1017: Level: beginner
1019: Notes:
1020: If the `Vec` was obtained from a `DM` with `DMCreateGlobalVector()`, then the range values are determined by the specific `DM`.
1022: If the `Vec` was created directly the range values are determined by the local size passed to `VecSetSizes()` or `VecCreateMPI()`.
1023: If `PETSC_DECIDE` was passed as the local size, then the vector uses default values for the range using `PetscSplitOwnership()`.
1025: The high argument is one more than the last element stored locally.
1027: For certain `DM`, such as `DMDA`, it is better to use `DM` specific routines, such as `DMDAGetGhostCorners()`, to determine
1028: the local values in the vector.
1030: .seealso: [](ch_vectors), `Vec`, `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `VecGetOwnershipRanges()`, `PetscSplitOwnership()`,
1031: `VecSetSizes()`, `VecCreateMPI()`, `PetscLayout`, `DMDAGetGhostCorners()`, `DM`
1032: @*/
1033: PetscErrorCode VecGetOwnershipRange(Vec x, PetscInt *low, PetscInt *high)
1034: {
1035: PetscFunctionBegin;
1038: if (low) PetscAssertPointer(low, 2);
1039: if (high) PetscAssertPointer(high, 3);
1040: if (low) *low = x->map->rstart;
1041: if (high) *high = x->map->rend;
1042: PetscFunctionReturn(PETSC_SUCCESS);
1043: }
1045: /*@C
1046: VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
1047: The vector is laid out with the
1048: first `n1` elements on the first processor, next `n2` elements on the
1049: second, etc. For certain parallel layouts this range may not be
1050: well defined.
1052: Not Collective
1054: Input Parameter:
1055: . x - the vector
1057: Output Parameter:
1058: . ranges - array of length `size` + 1 with the start and end+1 for each process
1060: Level: beginner
1062: Notes:
1063: If the `Vec` was obtained from a `DM` with `DMCreateGlobalVector()`, then the range values are determined by the specific `DM`.
1065: If the `Vec` was created directly the range values are determined by the local size passed to `VecSetSizes()` or `VecCreateMPI()`.
1066: If `PETSC_DECIDE` was passed as the local size, then the vector uses default values for the range using `PetscSplitOwnership()`.
1068: The high argument is one more than the last element stored locally.
1070: For certain `DM`, such as `DMDA`, it is better to use `DM` specific routines, such as `DMDAGetGhostCorners()`, to determine
1071: the local values in the vector.
1073: The high argument is one more than the last element stored locally.
1075: If `ranges` are used after all vectors that share the ranges has been destroyed, then the program will crash accessing `ranges`.
1077: Fortran Note:
1078: The argument `ranges` must be declared as
1079: .vb
1080: PetscInt, pointer :: ranges(:)
1081: .ve
1082: and you have to return it with a call to `VecRestoreOwnershipRanges()` when no longer needed
1084: .seealso: [](ch_vectors), `Vec`, `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `VecGetOwnershipRange()`, `PetscSplitOwnership()`,
1085: `VecSetSizes()`, `VecCreateMPI()`, `PetscLayout`, `DMDAGetGhostCorners()`, `DM`
1086: @*/
1087: PetscErrorCode VecGetOwnershipRanges(Vec x, const PetscInt *ranges[])
1088: {
1089: PetscFunctionBegin;
1092: PetscCall(PetscLayoutGetRanges(x->map, ranges));
1093: PetscFunctionReturn(PETSC_SUCCESS);
1094: }
1096: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
1097: /*@
1098: VecSetOption - Sets an option for controlling a vector's behavior.
1100: Collective
1102: Input Parameters:
1103: + x - the vector
1104: . op - the option
1105: - flag - turn the option on or off
1107: Supported Options:
1108: + `VEC_IGNORE_OFF_PROC_ENTRIES` - which causes `VecSetValues()` to ignore
1109: entries destined to be stored on a separate processor. This can be used
1110: to eliminate the global reduction in the `VecAssemblyBegin()` if you know
1111: that you have only used `VecSetValues()` to set local elements
1112: . `VEC_IGNORE_NEGATIVE_INDICES` - which means you can pass negative indices
1113: in ix in calls to `VecSetValues()` or `VecGetValues()`. These rows are simply
1114: ignored.
1115: - `VEC_SUBSET_OFF_PROC_ENTRIES` - which causes `VecAssemblyBegin()` to assume that the off-process
1116: entries will always be a subset (possibly equal) of the off-process entries set on the
1117: first assembly which had a true `VEC_SUBSET_OFF_PROC_ENTRIES` and the vector has not
1118: changed this flag afterwards. If this assembly is not such first assembly, then this
1119: assembly can reuse the communication pattern setup in that first assembly, thus avoiding
1120: a global reduction. Subsequent assemblies setting off-process values should use the same
1121: InsertMode as the first assembly.
1123: Level: intermediate
1125: Developer Notes:
1126: The `InsertMode` restriction could be removed by packing the stash messages out of place.
1128: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`
1129: @*/
1130: PetscErrorCode VecSetOption(Vec x, VecOption op, PetscBool flag)
1131: {
1132: PetscFunctionBegin;
1135: PetscTryTypeMethod(x, setoption, op, flag);
1136: PetscFunctionReturn(PETSC_SUCCESS);
1137: }
1139: /* Default routines for obtaining and releasing; */
1140: /* may be used by any implementation */
1141: PetscErrorCode VecDuplicateVecs_Default(Vec w, PetscInt m, Vec *V[])
1142: {
1143: PetscFunctionBegin;
1144: PetscCheck(m > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "m must be > 0: m = %" PetscInt_FMT, m);
1145: PetscCall(PetscMalloc1(m, V));
1146: for (PetscInt i = 0; i < m; i++) PetscCall(VecDuplicate(w, *V + i));
1147: PetscFunctionReturn(PETSC_SUCCESS);
1148: }
1150: PetscErrorCode VecDestroyVecs_Default(PetscInt m, Vec v[])
1151: {
1152: PetscInt i;
1154: PetscFunctionBegin;
1155: PetscAssertPointer(v, 2);
1156: for (i = 0; i < m; i++) PetscCall(VecDestroy(&v[i]));
1157: PetscCall(PetscFree(v));
1158: PetscFunctionReturn(PETSC_SUCCESS);
1159: }
1161: /*@
1162: VecResetArray - Resets a vector to use its default memory. Call this
1163: after the use of `VecPlaceArray()`.
1165: Not Collective
1167: Input Parameter:
1168: . vec - the vector
1170: Level: developer
1172: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecPlaceArray()`
1173: @*/
1174: PetscErrorCode VecResetArray(Vec vec)
1175: {
1176: PetscFunctionBegin;
1179: PetscUseTypeMethod(vec, resetarray);
1180: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
1181: PetscFunctionReturn(PETSC_SUCCESS);
1182: }
1184: /*@
1185: VecLoad - Loads a vector that has been stored in binary or HDF5 format
1186: with `VecView()`.
1188: Collective
1190: Input Parameters:
1191: + vec - the newly loaded vector, this needs to have been created with `VecCreate()` or
1192: some related function before the call to `VecLoad()`.
1193: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()` or
1194: HDF5 file viewer, obtained from `PetscViewerHDF5Open()`
1196: Level: intermediate
1198: Notes:
1199: Defaults to the standard `VECSEQ` or `VECMPI`, if you want some other type of `Vec` call `VecSetFromOptions()`
1200: before calling this.
1202: The input file must contain the full global vector, as
1203: written by the routine `VecView()`.
1205: If the type or size of `vec` is not set before a call to `VecLoad()`, PETSc
1206: sets the type and the local and global sizes based on the vector it is reading in. If type and/or
1207: sizes are already set, then the same are used.
1209: If using the binary viewer and the blocksize of the vector is greater than one then you must provide a unique prefix to
1210: the vector with `PetscObjectSetOptionsPrefix`((`PetscObject`)vec,"uniqueprefix"); BEFORE calling `VecView()` on the
1211: vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
1212: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
1213: filename. If you copy the binary file, make sure you copy the associated .info file with it.
1215: If using HDF5, you must assign the `Vec` the same name as was used in the Vec
1216: that was stored in the file using `PetscObjectSetName()`. Otherwise you will
1217: get the error message: "Cannot H5DOpen2() with `Vec` name NAMEOFOBJECT".
1219: If the HDF5 file contains a two dimensional array the first dimension is treated as the block size
1220: in loading the vector. Hence, for example, using MATLAB notation h5create('vector.dat','/Test_Vec',[27 1]);
1221: will load a vector of size 27 and block size 27 thus resulting in all 27 entries being on the first process of
1222: vectors communicator and the rest of the processes having zero entries
1224: Notes for advanced users when using the binary viewer:
1225: Most users should not need to know the details of the binary storage
1226: format, since `VecLoad()` and `VecView()` completely hide these details.
1227: But for anyone who's interested, the standard binary vector storage
1228: format is
1229: .vb
1230: PetscInt VEC_FILE_CLASSID
1231: PetscInt number of rows
1232: PetscScalar *values of all entries
1233: .ve
1235: In addition, PETSc automatically uses byte swapping to work on all machines; the files
1236: are written ALWAYS using big-endian ordering. On small-endian machines the numbers
1237: are converted to the small-endian format when they are read in from the file.
1238: See PetscBinaryRead() and PetscBinaryWrite() to see how this may be done.
1240: .seealso: [](ch_vectors), `Vec`, `PetscViewerBinaryOpen()`, `VecView()`, `MatLoad()`
1241: @*/
1242: PetscErrorCode VecLoad(Vec vec, PetscViewer viewer)
1243: {
1244: PetscViewerFormat format;
1246: PetscFunctionBegin;
1249: PetscCheckSameComm(vec, 1, viewer, 2);
1251: PetscCall(VecSetErrorIfLocked(vec, 1));
1252: if (!((PetscObject)vec)->type_name && !vec->ops->create) PetscCall(VecSetType(vec, VECSTANDARD));
1253: PetscCall(PetscLogEventBegin(VEC_Load, viewer, 0, 0, 0));
1254: PetscCall(PetscViewerGetFormat(viewer, &format));
1255: if (format == PETSC_VIEWER_NATIVE && vec->ops->loadnative) {
1256: PetscUseTypeMethod(vec, loadnative, viewer);
1257: } else {
1258: PetscUseTypeMethod(vec, load, viewer);
1259: }
1260: PetscCall(PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0));
1261: PetscFunctionReturn(PETSC_SUCCESS);
1262: }
1264: /*@
1265: VecReciprocal - Replaces each component of a vector by its reciprocal.
1267: Logically Collective
1269: Input Parameter:
1270: . vec - the vector
1272: Output Parameter:
1273: . vec - the vector reciprocal
1275: Level: intermediate
1277: Note:
1278: Vector entries with value 0.0 are not changed
1280: .seealso: [](ch_vectors), `Vec`, `VecLog()`, `VecExp()`, `VecSqrtAbs()`
1281: @*/
1282: PetscErrorCode VecReciprocal(Vec vec)
1283: {
1284: PetscFunctionBegin;
1285: PetscCall(VecReciprocalAsync_Private(vec, NULL));
1286: PetscFunctionReturn(PETSC_SUCCESS);
1287: }
1289: /*@C
1290: VecSetOperation - Allows the user to override a particular vector operation.
1292: Logically Collective; No Fortran Support
1294: Input Parameters:
1295: + vec - The vector to modify
1296: . op - The name of the operation
1297: - f - The function that provides the operation.
1299: Level: advanced
1301: Example Usage:
1302: .vb
1303: // some new VecView() implementation, must have the same signature as the function it seeks
1304: // to replace
1305: PetscErrorCode UserVecView(Vec x, PetscViewer viewer)
1306: {
1307: PetscFunctionBeginUser;
1308: // ...
1309: PetscFunctionReturn(PETSC_SUCCESS);
1310: }
1312: // Create a VECMPI which has a pre-defined VecView() implementation
1313: VecCreateMPI(comm, n, N, &x);
1314: // Calls the VECMPI implementation for VecView()
1315: VecView(x, viewer);
1317: VecSetOperation(x, VECOP_VIEW, (PetscErrorCodeFn *)UserVecView);
1318: // Now calls UserVecView()
1319: VecView(x, viewer);
1320: .ve
1322: Notes:
1323: `f` may be `NULL` to remove the operation from `vec`. Depending on the operation this may be
1324: allowed, however some always expect a valid function. In these cases an error will be raised
1325: when calling the interface routine in question.
1327: See `VecOperation` for an up-to-date list of override-able operations. The operations listed
1328: there have the form `VECOP_<OPERATION>`, where `<OPERATION>` is the suffix (in all capital
1329: letters) of the public interface routine (e.g., `VecView()` -> `VECOP_VIEW`).
1331: Overriding a particular `Vec`'s operation has no affect on any other `Vec`s past, present,
1332: or future. The user should also note that overriding a method is "destructive"; the previous
1333: method is not retained in any way.
1335: Each function MUST return `PETSC_SUCCESS` on success and
1336: nonzero on failure.
1338: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecGetOperation()`, `MatSetOperation()`, `MatShellSetOperation()`
1339: @*/
1340: PetscErrorCode VecSetOperation(Vec vec, VecOperation op, PetscErrorCodeFn *f)
1341: {
1342: PetscFunctionBegin;
1344: if (op == VECOP_VIEW && !vec->ops->viewnative) {
1345: vec->ops->viewnative = vec->ops->view;
1346: } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
1347: vec->ops->loadnative = vec->ops->load;
1348: }
1349: ((PetscErrorCodeFn **)vec->ops)[(int)op] = f;
1350: PetscFunctionReturn(PETSC_SUCCESS);
1351: }
1353: /*@
1354: VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1355: used during the assembly process to store values that belong to
1356: other processors.
1358: Not Collective, different processes can have different size stashes
1360: Input Parameters:
1361: + vec - the vector
1362: . size - the initial size of the stash.
1363: - bsize - the initial size of the block-stash(if used).
1365: Options Database Keys:
1366: + -vecstash_initial_size size or size0,size1,...,sizep-1 - set initial size
1367: - -vecstash_block_initial_size bsize or bsize0,bsize1,...,bsizep-1 - set initial block size
1369: Level: intermediate
1371: Notes:
1372: The block-stash is used for values set with `VecSetValuesBlocked()` while
1373: the stash is used for values set with `VecSetValues()`
1375: Run with the option -info and look for output of the form
1376: VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1377: to determine the appropriate value, MM, to use for size and
1378: VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1379: to determine the value, BMM to use for bsize
1381: PETSc attempts to smartly manage the stash size so there is little likelihood setting a
1382: a specific value here will affect performance
1384: .seealso: [](ch_vectors), `Vec`, `VecSetBlockSize()`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecStashView()`
1385: @*/
1386: PetscErrorCode VecStashSetInitialSize(Vec vec, PetscInt size, PetscInt bsize)
1387: {
1388: PetscFunctionBegin;
1390: PetscCall(VecStashSetInitialSize_Private(&vec->stash, size));
1391: PetscCall(VecStashSetInitialSize_Private(&vec->bstash, bsize));
1392: PetscFunctionReturn(PETSC_SUCCESS);
1393: }
1395: /*@
1396: VecSetRandom - Sets all components of a vector to random numbers.
1398: Logically Collective
1400: Input Parameters:
1401: + x - the vector
1402: - rctx - the random number context, formed by `PetscRandomCreate()`, or use `NULL` and it will create one internally.
1404: Output Parameter:
1405: . x - the vector
1407: Example of Usage:
1408: .vb
1409: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1410: VecSetRandom(x,rctx);
1411: PetscRandomDestroy(&rctx);
1412: .ve
1414: Level: intermediate
1416: .seealso: [](ch_vectors), `Vec`, `VecSet()`, `VecSetValues()`, `PetscRandomCreate()`, `PetscRandomDestroy()`
1417: @*/
1418: PetscErrorCode VecSetRandom(Vec x, PetscRandom rctx)
1419: {
1420: PetscRandom randObj = NULL;
1422: PetscFunctionBegin;
1426: VecCheckAssembled(x);
1427: PetscCall(VecSetErrorIfLocked(x, 1));
1429: if (!rctx) {
1430: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)x), &randObj));
1431: PetscCall(PetscRandomSetType(randObj, x->defaultrandtype));
1432: PetscCall(PetscRandomSetFromOptions(randObj));
1433: rctx = randObj;
1434: }
1436: PetscCall(PetscLogEventBegin(VEC_SetRandom, x, rctx, 0, 0));
1437: PetscUseTypeMethod(x, setrandom, rctx);
1438: PetscCall(PetscLogEventEnd(VEC_SetRandom, x, rctx, 0, 0));
1440: PetscCall(PetscRandomDestroy(&randObj));
1441: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1442: PetscFunctionReturn(PETSC_SUCCESS);
1443: }
1445: /*@
1446: VecSetRandomGaussian - Fills a vector with Gaussian random values of the given mean and standard deviation.
1448: Collective
1450: Input Parameters:
1451: + v - the vector to fill
1452: . rng - PETSc random number generator
1453: . mean - desired mean of the Gaussian samples
1454: - std_dev - desired standard deviation
1456: Level: advanced
1458: Note:
1459: For complex builds where `PetscScalar` is complex the imaginary part of all the vector entries is zero
1461: Developer Note:
1462: Uses the Box-Muller transform to generate normally distributed random numbers
1463: from uniform random numbers. Handles edge cases where uniform random values
1464: approach 0 or 1.
1466: .seealso: [](ch_vectors), [](ch_da), `PetscDA`, `PetscRandom`, `PetscRandomSetInterval()`, `VecSetRandom()`
1467: @*/
1468: PetscErrorCode VecSetRandomGaussian(Vec v, PetscRandom rng, PetscReal mean, PetscReal std_dev)
1469: {
1470: PetscInt n, i;
1471: PetscScalar *array;
1472: PetscReal u1, u2;
1473: PetscReal gauss_sample1, gauss_sample2, magnitude, theta;
1474: const PetscReal min_uniform = PETSC_MACHINE_EPSILON;
1475: const PetscInt max_retry_count = 100;
1477: PetscFunctionBegin;
1482: PetscCheck(!PetscIsInfOrNanReal(mean), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Mean must be a finite real number");
1483: PetscCheck(std_dev >= 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Standard deviation must be non-negative, got %g", (double)std_dev);
1484: PetscCheck(!PetscIsInfOrNanReal(std_dev), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Standard deviation must be a finite real number");
1486: PetscCall(VecGetLocalSize(v, &n));
1487: if (n == 0) PetscFunctionReturn(PETSC_SUCCESS);
1489: if (std_dev == 0.0) {
1490: PetscCall(VecSet(v, mean));
1491: PetscFunctionReturn(PETSC_SUCCESS);
1492: }
1494: PetscCall(VecGetArrayWrite(v, &array));
1496: /*
1497: Generate Gaussian-distributed random values using the Box-Muller transform.
1498: This transform converts pairs of uniform random variables U1, U2 ~ Uniform(0,1)
1499: into pairs of independent standard normal variables Z0, Z1 ~ N(0,1):
1500: Z0 = sqrt(-2 * ln(U1)) * cos(2pi * U2)
1501: Z1 = sqrt(-2 * ln(U1)) * sin(2pi * U2)
1502: Then scale and shift to get desired mean and standard deviation.
1503: */
1504: for (i = 0; i < n; i += 2) {
1505: PetscInt retry_count = 0;
1507: /*
1508: Generate U1 and ensure it's not too close to 0 to avoid log(0) singularity.
1509: Add retry limit to prevent infinite loops in case of RNG failure.
1510: */
1511: do {
1512: PetscCall(PetscRandomGetValueReal(rng, &u1));
1513: retry_count++;
1514: 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);
1515: } while (u1 < min_uniform);
1517: PetscCall(PetscRandomGetValueReal(rng, &u2));
1519: /*
1520: Apply Box-Muller transform:
1521: - magnitude: sqrt(-2 * ln(U1)) represents the radial distance from origin
1522: - theta: 2pi * U2 represents the angle uniformly distributed on [0, 2pi]
1523: - Converting from polar to Cartesian coordinates yields two independent samples
1524: */
1525: magnitude = PetscSqrtReal(-2.0 * PetscLogReal(u1));
1526: theta = 2.0 * PETSC_PI * u2;
1527: gauss_sample1 = magnitude * PetscCosReal(theta);
1528: gauss_sample2 = magnitude * PetscSinReal(theta);
1530: /* Scale and shift to achieve desired mean and standard deviation */
1531: array[i] = mean + std_dev * gauss_sample1;
1532: if (i + 1 < n) array[i + 1] = mean + std_dev * gauss_sample2;
1533: }
1535: PetscCall(VecRestoreArrayWrite(v, &array));
1536: PetscFunctionReturn(PETSC_SUCCESS);
1537: }
1539: /*@
1540: VecZeroEntries - puts a `0.0` in each element of a vector
1542: Logically Collective
1544: Input Parameter:
1545: . vec - The vector
1547: Level: beginner
1549: Note:
1550: If the norm of the vector is known to be zero then this skips the unneeded zeroing process
1552: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecSetOptionsPrefix()`, `VecSet()`, `VecSetValues()`
1553: @*/
1554: PetscErrorCode VecZeroEntries(Vec vec)
1555: {
1556: PetscFunctionBegin;
1557: PetscCall(VecSet(vec, 0));
1558: PetscFunctionReturn(PETSC_SUCCESS);
1559: }
1561: /*
1562: VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
1563: processor and a PETSc MPI vector on more than one processor.
1565: Collective
1567: Input Parameter:
1568: . vec - The vector
1570: Level: intermediate
1572: .seealso: [](ch_vectors), `Vec`, `VecSetFromOptions()`, `VecSetType()`
1573: */
1574: static PetscErrorCode VecSetTypeFromOptions_Private(Vec vec, PetscOptionItems PetscOptionsObject)
1575: {
1576: PetscBool opt;
1577: VecType defaultType;
1578: char typeName[256];
1579: PetscMPIInt size;
1581: PetscFunctionBegin;
1582: if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1583: else {
1584: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
1585: if (size > 1) defaultType = VECMPI;
1586: else defaultType = VECSEQ;
1587: }
1589: PetscCall(VecRegisterAll());
1590: PetscCall(PetscOptionsFList("-vec_type", "Vector type", "VecSetType", VecList, defaultType, typeName, 256, &opt));
1591: if (opt) PetscCall(VecSetType(vec, typeName));
1592: else PetscCall(VecSetType(vec, defaultType));
1593: PetscFunctionReturn(PETSC_SUCCESS);
1594: }
1596: /*@
1597: VecSetFromOptions - Configures the vector from the options database.
1599: Collective
1601: Input Parameter:
1602: . vec - The vector
1604: Level: beginner
1606: Notes:
1607: To see all options, run your program with the -help option.
1609: Must be called after `VecCreate()` but before the vector is used.
1611: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecSetOptionsPrefix()`
1612: @*/
1613: PetscErrorCode VecSetFromOptions(Vec vec)
1614: {
1615: PetscBool flg;
1616: PetscInt bind_below = 0;
1618: PetscFunctionBegin;
1621: PetscObjectOptionsBegin((PetscObject)vec);
1622: /* Handle vector type options */
1623: PetscCall(VecSetTypeFromOptions_Private(vec, PetscOptionsObject));
1625: /* Handle specific vector options */
1626: PetscTryTypeMethod(vec, setfromoptions, PetscOptionsObject);
1628: /* Bind to CPU if below a user-specified size threshold.
1629: * This perhaps belongs in the options for the GPU Vec types, but VecBindToCPU() does nothing when called on non-GPU types,
1630: * and putting it here makes is more maintainable than duplicating this for all. */
1631: 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));
1632: if (flg && vec->map->n < bind_below) PetscCall(VecBindToCPU(vec, PETSC_TRUE));
1634: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1635: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)vec, PetscOptionsObject));
1636: PetscOptionsEnd();
1637: PetscFunctionReturn(PETSC_SUCCESS);
1638: }
1640: /*@
1641: VecSetSizes - Sets the local and global sizes, and checks to determine compatibility of the sizes
1643: Collective
1645: Input Parameters:
1646: + v - the vector
1647: . n - the local size (or `PETSC_DECIDE` to have it set)
1648: - N - the global size (or `PETSC_DETERMINE` to have it set)
1650: Level: intermediate
1652: Notes:
1653: `N` cannot be `PETSC_DETERMINE` if `n` is `PETSC_DECIDE`
1655: If one processor calls this with `N` of `PETSC_DETERMINE` then all processors must, otherwise the program will hang.
1657: If `n` is not `PETSC_DECIDE`, then the value determines the `PetscLayout` of the vector and the ranges returned by
1658: `VecGetOwnershipRange()` and `VecGetOwnershipRanges()`
1660: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecCreateSeq()`, `VecCreateMPI()`, `VecGetSize()`, `PetscSplitOwnership()`, `PetscLayout`,
1661: `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`, `MatSetSizes()`
1662: @*/
1663: PetscErrorCode VecSetSizes(Vec v, PetscInt n, PetscInt N)
1664: {
1665: PetscFunctionBegin;
1667: if (N >= 0) {
1669: PetscCheck(n <= N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local size %" PetscInt_FMT " cannot be larger than global size %" PetscInt_FMT, n, N);
1670: }
1671: 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,
1672: v->map->n, v->map->N);
1673: v->map->n = n;
1674: v->map->N = N;
1675: PetscTryTypeMethod(v, create);
1676: v->ops->create = NULL;
1677: PetscFunctionReturn(PETSC_SUCCESS);
1678: }
1680: /*@
1681: VecSetBlockSize - Sets the block size for future calls to `VecSetValuesBlocked()`
1682: and `VecSetValuesBlockedLocal()`.
1684: Logically Collective
1686: Input Parameters:
1687: + v - the vector
1688: - bs - the blocksize
1690: Level: advanced
1692: Note:
1693: All vectors obtained by `VecDuplicate()` inherit the same blocksize.
1695: Vectors obtained with `DMCreateGlobalVector()` and `DMCreateLocalVector()` generally already have a blocksize set based on the state of the `DM`
1697: .seealso: [](ch_vectors), `Vec`, `VecSetValuesBlocked()`, `VecSetLocalToGlobalMapping()`, `VecGetBlockSize()`
1698: @*/
1699: PetscErrorCode VecSetBlockSize(Vec v, PetscInt bs)
1700: {
1701: PetscFunctionBegin;
1704: PetscCall(PetscLayoutSetBlockSize(v->map, bs));
1705: v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1706: PetscFunctionReturn(PETSC_SUCCESS);
1707: }
1709: /*@
1710: VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for `VecSetValuesBlocked()`
1711: and `VecSetValuesBlockedLocal()`.
1713: Not Collective
1715: Input Parameter:
1716: . v - the vector
1718: Output Parameter:
1719: . bs - the blocksize
1721: Level: advanced
1723: Note:
1724: All vectors obtained by `VecDuplicate()` inherit the same blocksize.
1726: .seealso: [](ch_vectors), `Vec`, `VecSetValuesBlocked()`, `VecSetLocalToGlobalMapping()`, `VecSetBlockSize()`
1727: @*/
1728: PetscErrorCode VecGetBlockSize(Vec v, PetscInt *bs)
1729: {
1730: PetscFunctionBegin;
1732: PetscAssertPointer(bs, 2);
1733: PetscCall(PetscLayoutGetBlockSize(v->map, bs));
1734: PetscFunctionReturn(PETSC_SUCCESS);
1735: }
1737: /*@
1738: VecSetOptionsPrefix - Sets the prefix used for searching for all
1739: `Vec` options in the database.
1741: Logically Collective
1743: Input Parameters:
1744: + v - the `Vec` context
1745: - prefix - the prefix to prepend to all option names
1747: Level: advanced
1749: Note:
1750: A hyphen (-) must NOT be given at the beginning of the prefix name.
1751: The first character of all runtime options is AUTOMATICALLY the hyphen.
1753: .seealso: [](ch_vectors), `Vec`, `VecSetFromOptions()`
1754: @*/
1755: PetscErrorCode VecSetOptionsPrefix(Vec v, const char prefix[])
1756: {
1757: PetscFunctionBegin;
1759: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)v, prefix));
1760: PetscFunctionReturn(PETSC_SUCCESS);
1761: }
1763: /*@
1764: VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1765: `Vec` options in the database.
1767: Logically Collective
1769: Input Parameters:
1770: + v - the `Vec` context
1771: - prefix - the prefix to prepend to all option names
1773: Level: advanced
1775: Note:
1776: A hyphen (-) must NOT be given at the beginning of the prefix name.
1777: The first character of all runtime options is AUTOMATICALLY the hyphen.
1779: .seealso: [](ch_vectors), `Vec`, `VecGetOptionsPrefix()`
1780: @*/
1781: PetscErrorCode VecAppendOptionsPrefix(Vec v, const char prefix[])
1782: {
1783: PetscFunctionBegin;
1785: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)v, prefix));
1786: PetscFunctionReturn(PETSC_SUCCESS);
1787: }
1789: /*@
1790: VecGetOptionsPrefix - Sets the prefix used for searching for all
1791: Vec options in the database.
1793: Not Collective
1795: Input Parameter:
1796: . v - the `Vec` context
1798: Output Parameter:
1799: . prefix - pointer to the prefix string used
1801: Level: advanced
1803: .seealso: [](ch_vectors), `Vec`, `VecAppendOptionsPrefix()`
1804: @*/
1805: PetscErrorCode VecGetOptionsPrefix(Vec v, const char *prefix[])
1806: {
1807: PetscFunctionBegin;
1809: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)v, prefix));
1810: PetscFunctionReturn(PETSC_SUCCESS);
1811: }
1813: /*@C
1814: VecGetState - Gets the state of a `Vec`.
1816: Not Collective
1818: Input Parameter:
1819: . v - the `Vec` context
1821: Output Parameter:
1822: . state - the object state
1824: Level: advanced
1826: Note:
1827: Object state is an integer which gets increased every time
1828: the object is changed. By saving and later querying the object state
1829: one can determine whether information about the object is still current.
1831: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `PetscObjectStateGet()`
1832: @*/
1833: PetscErrorCode VecGetState(Vec v, PetscObjectState *state)
1834: {
1835: PetscFunctionBegin;
1837: PetscAssertPointer(state, 2);
1838: PetscCall(PetscObjectStateGet((PetscObject)v, state));
1839: PetscFunctionReturn(PETSC_SUCCESS);
1840: }
1842: /*@
1843: VecSetUp - Sets up the internal vector data structures for the later use.
1845: Collective
1847: Input Parameter:
1848: . v - the `Vec` context
1850: Level: advanced
1852: Notes:
1853: For basic use of the `Vec` classes the user need not explicitly call
1854: `VecSetUp()`, since these actions will happen automatically.
1856: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecDestroy()`
1857: @*/
1858: PetscErrorCode VecSetUp(Vec v)
1859: {
1860: PetscMPIInt size;
1862: PetscFunctionBegin;
1864: PetscCheck(v->map->n >= 0 || v->map->N >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Sizes not set");
1865: if (!((PetscObject)v)->type_name) {
1866: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1867: if (size == 1) PetscCall(VecSetType(v, VECSEQ));
1868: else PetscCall(VecSetType(v, VECMPI));
1869: }
1870: PetscFunctionReturn(PETSC_SUCCESS);
1871: }
1873: /*
1874: These currently expose the PetscScalar/PetscReal in updating the
1875: cached norm. If we push those down into the implementation these
1876: will become independent of PetscScalar/PetscReal
1877: */
1879: PetscErrorCode VecCopyAsync_Private(Vec x, Vec y, PetscDeviceContext dctx)
1880: {
1881: PetscBool flgs[4];
1882: PetscReal norms[4] = {0.0, 0.0, 0.0, 0.0};
1884: PetscFunctionBegin;
1889: if (x == y) PetscFunctionReturn(PETSC_SUCCESS);
1890: VecCheckSameLocalSize(x, 1, y, 2);
1891: VecCheckAssembled(x);
1892: PetscCall(VecSetErrorIfLocked(y, 2));
1894: #if !defined(PETSC_USE_MIXED_PRECISION)
1895: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
1896: #endif
1898: PetscCall(PetscLogEventBegin(VEC_Copy, x, y, 0, 0));
1899: #if defined(PETSC_USE_MIXED_PRECISION)
1900: extern PetscErrorCode VecGetArray(Vec, double **);
1901: extern PetscErrorCode VecRestoreArray(Vec, double **);
1902: extern PetscErrorCode VecGetArray(Vec, float **);
1903: extern PetscErrorCode VecRestoreArray(Vec, float **);
1904: extern PetscErrorCode VecGetArrayRead(Vec, const double **);
1905: extern PetscErrorCode VecRestoreArrayRead(Vec, const double **);
1906: extern PetscErrorCode VecGetArrayRead(Vec, const float **);
1907: extern PetscErrorCode VecRestoreArrayRead(Vec, const float **);
1908: if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1909: PetscInt i, n;
1910: const float *xx;
1911: double *yy;
1912: PetscCall(VecGetArrayRead(x, &xx));
1913: PetscCall(VecGetArray(y, &yy));
1914: PetscCall(VecGetLocalSize(x, &n));
1915: for (i = 0; i < n; i++) yy[i] = xx[i];
1916: PetscCall(VecRestoreArrayRead(x, &xx));
1917: PetscCall(VecRestoreArray(y, &yy));
1918: } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1919: PetscInt i, n;
1920: float *yy;
1921: const double *xx;
1922: PetscCall(VecGetArrayRead(x, &xx));
1923: PetscCall(VecGetArray(y, &yy));
1924: PetscCall(VecGetLocalSize(x, &n));
1925: for (i = 0; i < n; i++) yy[i] = (float)xx[i];
1926: PetscCall(VecRestoreArrayRead(x, &xx));
1927: PetscCall(VecRestoreArray(y, &yy));
1928: } else PetscUseTypeMethod(x, copy, y);
1929: #else
1930: VecMethodDispatch(x, dctx, VecAsyncFnName(Copy), copy, (Vec, Vec, PetscDeviceContext), y);
1931: #endif
1933: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1934: #if !defined(PETSC_USE_MIXED_PRECISION)
1935: for (PetscInt i = 0; i < 4; i++) {
1936: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)y, NormIds[i], norms[i]));
1937: }
1938: #endif
1940: PetscCall(PetscLogEventEnd(VEC_Copy, x, y, 0, 0));
1941: PetscFunctionReturn(PETSC_SUCCESS);
1942: }
1944: /*@
1945: VecCopy - Copies a vector `y = x`
1947: Logically Collective
1949: Input Parameter:
1950: . x - the vector
1952: Output Parameter:
1953: . y - the copy
1955: Level: beginner
1957: Note:
1958: For default parallel PETSc vectors, both `x` and `y` must be distributed in
1959: the same manner; local copies are done.
1961: Developer Notes:
1962: `PetscCheckSameTypeAndComm`(x,1,y,2) is not used on these vectors because we allow one
1963: of the vectors to be sequential and one to be parallel so long as both have the same
1964: local sizes. This is used in some internal functions in PETSc.
1966: .seealso: [](ch_vectors), `Vec`, `VecDuplicate()`
1967: @*/
1968: PetscErrorCode VecCopy(Vec x, Vec y)
1969: {
1970: PetscFunctionBegin;
1971: PetscCall(VecCopyAsync_Private(x, y, NULL));
1972: PetscFunctionReturn(PETSC_SUCCESS);
1973: }
1975: PetscErrorCode VecSwapAsync_Private(Vec x, Vec y, PetscDeviceContext dctx)
1976: {
1977: PetscReal normxs[4], normys[4];
1978: PetscBool flgxs[4], flgys[4];
1980: PetscFunctionBegin;
1985: PetscCheckSameTypeAndComm(x, 1, y, 2);
1986: VecCheckSameSize(x, 1, y, 2);
1987: VecCheckAssembled(x);
1988: VecCheckAssembled(y);
1989: PetscCall(VecSetErrorIfLocked(x, 1));
1990: PetscCall(VecSetErrorIfLocked(y, 2));
1992: for (PetscInt i = 0; i < 4; i++) {
1993: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], normxs[i], flgxs[i]));
1994: PetscCall(PetscObjectComposedDataGetReal((PetscObject)y, NormIds[i], normys[i], flgys[i]));
1995: }
1997: PetscCall(PetscLogEventBegin(VEC_Swap, x, y, 0, 0));
1998: VecMethodDispatch(x, dctx, VecAsyncFnName(Swap), swap, (Vec, Vec, PetscDeviceContext), y);
1999: PetscCall(PetscLogEventEnd(VEC_Swap, x, y, 0, 0));
2001: PetscCall(PetscObjectStateIncrease((PetscObject)x));
2002: PetscCall(PetscObjectStateIncrease((PetscObject)y));
2003: for (PetscInt i = 0; i < 4; i++) {
2004: if (flgxs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)y, NormIds[i], normxs[i]));
2005: if (flgys[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], normys[i]));
2006: }
2007: PetscFunctionReturn(PETSC_SUCCESS);
2008: }
2009: /*@
2010: VecSwap - Swaps the values between two vectors, `x` and `y`.
2012: Logically Collective
2014: Input Parameters:
2015: + x - the first vector
2016: - y - the second vector
2018: Level: advanced
2020: .seealso: [](ch_vectors), `Vec`, `VecSet()`
2021: @*/
2022: PetscErrorCode VecSwap(Vec x, Vec y)
2023: {
2024: PetscFunctionBegin;
2025: PetscCall(VecSwapAsync_Private(x, y, NULL));
2026: PetscFunctionReturn(PETSC_SUCCESS);
2027: }
2029: /*@
2030: VecStashViewFromOptions - Processes command line options to determine if/how a `VecStash` object is to be viewed.
2032: Collective
2034: Input Parameters:
2035: + obj - the `Vec` containing a stash
2036: . bobj - optional other object that provides the prefix
2037: - optionname - option to activate viewing
2039: Level: intermediate
2041: Developer Notes:
2042: This cannot use `PetscObjectViewFromOptions()` because it takes a `Vec` as an argument but does not use `VecView()`
2044: .seealso: [](ch_vectors), `Vec`, `VecStashSetInitialSize()`
2045: @*/
2046: PetscErrorCode VecStashViewFromOptions(Vec obj, PetscObject bobj, const char optionname[])
2047: {
2048: PetscViewer viewer;
2049: PetscBool flg;
2050: PetscViewerFormat format;
2051: char *prefix;
2053: PetscFunctionBegin;
2054: prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
2055: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)obj), ((PetscObject)obj)->options, prefix, optionname, &viewer, &format, &flg));
2056: if (flg) {
2057: PetscCall(PetscViewerPushFormat(viewer, format));
2058: PetscCall(VecStashView(obj, viewer));
2059: PetscCall(PetscViewerPopFormat(viewer));
2060: PetscCall(PetscViewerDestroy(&viewer));
2061: }
2062: PetscFunctionReturn(PETSC_SUCCESS);
2063: }
2065: /*@
2066: VecStashView - Prints the entries in the vector stash and block stash.
2068: Collective
2070: Input Parameters:
2071: + v - the vector
2072: - viewer - the viewer
2074: Level: advanced
2076: .seealso: [](ch_vectors), `Vec`, `VecSetBlockSize()`, `VecSetValues()`, `VecSetValuesBlocked()`
2077: @*/
2078: PetscErrorCode VecStashView(Vec v, PetscViewer viewer)
2079: {
2080: PetscMPIInt rank;
2081: PetscInt i, j;
2082: PetscBool match;
2083: VecStash *s;
2084: PetscScalar val;
2086: PetscFunctionBegin;
2089: PetscCheckSameComm(v, 1, viewer, 2);
2091: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &match));
2092: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_SUP, "Stash viewer only works with ASCII viewer not %s", ((PetscObject)v)->type_name);
2093: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
2094: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank));
2095: s = &v->bstash;
2097: /* print block stash */
2098: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2099: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Vector Block stash size %" PetscInt_FMT " block size %" PetscInt_FMT "\n", rank, s->n, s->bs));
2100: for (i = 0; i < s->n; i++) {
2101: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " ", rank, s->idx[i]));
2102: for (j = 0; j < s->bs; j++) {
2103: val = s->array[i * s->bs + j];
2104: #if defined(PETSC_USE_COMPLEX)
2105: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "(%18.16e %18.16e) ", (double)PetscRealPart(val), (double)PetscImaginaryPart(val)));
2106: #else
2107: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%18.16e ", (double)val));
2108: #endif
2109: }
2110: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
2111: }
2112: PetscCall(PetscViewerFlush(viewer));
2114: s = &v->stash;
2116: /* print basic stash */
2117: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Vector stash size %" PetscInt_FMT "\n", rank, s->n));
2118: for (i = 0; i < s->n; i++) {
2119: val = s->array[i];
2120: #if defined(PETSC_USE_COMPLEX)
2121: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " (%18.16e %18.16e) ", rank, s->idx[i], (double)PetscRealPart(val), (double)PetscImaginaryPart(val)));
2122: #else
2123: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " %18.16e\n", rank, s->idx[i], (double)val));
2124: #endif
2125: }
2126: PetscCall(PetscViewerFlush(viewer));
2127: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2128: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
2129: PetscFunctionReturn(PETSC_SUCCESS);
2130: }
2132: PetscErrorCode PetscOptionsGetVec(PetscOptions options, const char prefix[], const char key[], Vec v, PetscBool *set)
2133: {
2134: PetscInt i, N, rstart, rend;
2135: PetscScalar *xx;
2136: PetscReal *xreal;
2137: PetscBool iset;
2139: PetscFunctionBegin;
2140: PetscCall(VecGetOwnershipRange(v, &rstart, &rend));
2141: PetscCall(VecGetSize(v, &N));
2142: PetscCall(PetscCalloc1(N, &xreal));
2143: PetscCall(PetscOptionsGetRealArray(options, prefix, key, xreal, &N, &iset));
2144: if (iset) {
2145: PetscCall(VecGetArray(v, &xx));
2146: for (i = rstart; i < rend; i++) xx[i - rstart] = xreal[i];
2147: PetscCall(VecRestoreArray(v, &xx));
2148: }
2149: PetscCall(PetscFree(xreal));
2150: if (set) *set = iset;
2151: PetscFunctionReturn(PETSC_SUCCESS);
2152: }
2154: /*@
2155: VecGetLayout - get `PetscLayout` describing a vector layout
2157: Not Collective
2159: Input Parameter:
2160: . x - the vector
2162: Output Parameter:
2163: . map - the layout
2165: Level: developer
2167: Note:
2168: The layout determines what vector elements are contained on each MPI process
2170: .seealso: [](ch_vectors), `PetscLayout`, `Vec`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
2171: @*/
2172: PetscErrorCode VecGetLayout(Vec x, PetscLayout *map)
2173: {
2174: PetscFunctionBegin;
2176: PetscAssertPointer(map, 2);
2177: *map = x->map;
2178: PetscFunctionReturn(PETSC_SUCCESS);
2179: }
2181: /*@
2182: VecSetLayout - set `PetscLayout` describing vector layout
2184: Not Collective
2186: Input Parameters:
2187: + x - the vector
2188: - map - the layout
2190: Level: developer
2192: Note:
2193: It is normally only valid to replace the layout with a layout known to be equivalent.
2195: .seealso: [](ch_vectors), `Vec`, `PetscLayout`, `VecGetLayout()`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
2196: @*/
2197: PetscErrorCode VecSetLayout(Vec x, PetscLayout map)
2198: {
2199: PetscFunctionBegin;
2201: PetscCall(PetscLayoutReference(map, &x->map));
2202: PetscFunctionReturn(PETSC_SUCCESS);
2203: }
2205: /*@
2206: VecFlag - set infinity into the local part of the vector on any subset of MPI processes
2208: Logically Collective
2210: Input Parameters:
2211: + xin - the vector, can be `NULL` but only if on all processes
2212: - flg - indicates if this processes portion of the vector should be set to infinity
2214: Level: developer
2216: Note:
2217: This removes the values from the vector norm cache for all processes by calling `PetscObjectIncrease()`.
2219: This is used for any subset of MPI processes to indicate an failure in a solver, after the next use of `VecNorm()` if
2220: `KSPCheckNorm()` detects an infinity and at least one of the MPI processes has a not converged reason then the `KSP`
2221: object collectively is labeled as not converged.
2223: .seealso: [](ch_vectors), `Vec`, `PetscLayout`, `VecGetLayout()`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
2224: @*/
2225: PetscErrorCode VecFlag(Vec xin, PetscInt flg)
2226: {
2227: // MSVC gives "divide by zero" error at compile time - so declare as volatile to skip this check.
2228: volatile PetscReal one = 1.0, zero = 0.0;
2229: PetscScalar inf;
2231: PetscFunctionBegin;
2232: if (!xin) PetscFunctionReturn(PETSC_SUCCESS);
2234: PetscCall(PetscObjectStateIncrease((PetscObject)xin));
2235: if (flg) {
2236: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2237: inf = one / zero;
2238: PetscCall(PetscFPTrapPop());
2239: if (xin->ops->set) PetscUseTypeMethod(xin, set, inf);
2240: else {
2241: PetscInt n;
2242: PetscScalar *xx;
2244: PetscCall(VecGetLocalSize(xin, &n));
2245: PetscCall(VecGetArrayWrite(xin, &xx));
2246: for (PetscInt i = 0; i < n; ++i) xx[i] = inf;
2247: PetscCall(VecRestoreArrayWrite(xin, &xx));
2248: }
2249: }
2250: PetscFunctionReturn(PETSC_SUCCESS);
2251: }
2253: /*@
2254: VecSetInf - set infinity into the local part of the vector
2256: Not Collective
2258: Input Parameters:
2259: . xin - the vector
2261: Level: developer
2263: Note:
2264: Deprecated, see `VecFlag()`
2265: This is used for any subset of MPI processes to indicate an failure in a solver, after the next use of `VecNorm()` if
2266: `KSPCheckNorm()` detects an infinity and at least one of the MPI processes has a not converged reason then the `KSP`
2267: object collectively is labeled as not converged.
2269: This cannot be called if `xin` has a cached norm available
2271: .seealso: [](ch_vectors), `VecFlag()`, `Vec`, `PetscLayout`, `VecGetLayout()`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
2272: @*/
2273: PetscErrorCode VecSetInf(Vec xin)
2274: {
2275: // MSVC gives "divide by zero" error at compile time - so declare as volatile to skip this check.
2276: volatile PetscReal one = 1.0, zero = 0.0;
2277: PetscScalar inf;
2278: PetscBool flg;
2280: PetscFunctionBegin;
2281: PetscCall(VecNormAvailable(xin, NORM_2, &flg, NULL));
2282: PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot call VecSetInf() if the vector has a cached norm");
2283: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2284: inf = one / zero;
2285: PetscCall(PetscFPTrapPop());
2286: if (xin->ops->set) PetscUseTypeMethod(xin, set, inf);
2287: else {
2288: PetscInt n;
2289: PetscScalar *xx;
2291: PetscCall(VecGetLocalSize(xin, &n));
2292: PetscCall(VecGetArrayWrite(xin, &xx));
2293: for (PetscInt i = 0; i < n; ++i) xx[i] = inf;
2294: PetscCall(VecRestoreArrayWrite(xin, &xx));
2295: }
2296: PetscFunctionReturn(PETSC_SUCCESS);
2297: }
2299: /*@
2300: VecBindToCPU - marks a vector to temporarily stay on the CPU and perform computations on the CPU
2302: Logically collective
2304: Input Parameters:
2305: + v - the vector
2306: - flg - bind to the CPU if value of `PETSC_TRUE`
2308: Level: intermediate
2310: .seealso: [](ch_vectors), `Vec`, `VecBoundToCPU()`
2311: @*/
2312: PetscErrorCode VecBindToCPU(Vec v, PetscBool flg)
2313: {
2314: PetscFunctionBegin;
2317: #if defined(PETSC_HAVE_DEVICE)
2318: if (v->boundtocpu == flg) PetscFunctionReturn(PETSC_SUCCESS);
2319: v->boundtocpu = flg;
2320: PetscTryTypeMethod(v, bindtocpu, flg);
2321: #endif
2322: PetscFunctionReturn(PETSC_SUCCESS);
2323: }
2325: /*@
2326: VecBoundToCPU - query if a vector is bound to the CPU
2328: Not collective
2330: Input Parameter:
2331: . v - the vector
2333: Output Parameter:
2334: . flg - the logical flag
2336: Level: intermediate
2338: .seealso: [](ch_vectors), `Vec`, `VecBindToCPU()`
2339: @*/
2340: PetscErrorCode VecBoundToCPU(Vec v, PetscBool *flg)
2341: {
2342: PetscFunctionBegin;
2344: PetscAssertPointer(flg, 2);
2345: #if defined(PETSC_HAVE_DEVICE)
2346: *flg = v->boundtocpu;
2347: #else
2348: *flg = PETSC_TRUE;
2349: #endif
2350: PetscFunctionReturn(PETSC_SUCCESS);
2351: }
2353: /*@
2354: VecSetBindingPropagates - Sets whether the state of being bound to the CPU for a GPU vector type propagates to child and some other associated objects
2356: Input Parameters:
2357: + v - the vector
2358: - flg - flag indicating whether the boundtocpu flag should be propagated
2360: Level: developer
2362: Notes:
2363: 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.
2364: The created vectors will also have their bindingpropagates flag set to true.
2366: Developer Notes:
2367: If a `DMDA` has the `-dm_bind_below option` set to true, then vectors created by `DMCreateGlobalVector()` will have `VecSetBindingPropagates()` called on them to
2368: set their bindingpropagates flag to true.
2370: .seealso: [](ch_vectors), `Vec`, `MatSetBindingPropagates()`, `VecGetBindingPropagates()`
2371: @*/
2372: PetscErrorCode VecSetBindingPropagates(Vec v, PetscBool flg)
2373: {
2374: PetscFunctionBegin;
2376: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2377: v->bindingpropagates = flg;
2378: #endif
2379: PetscFunctionReturn(PETSC_SUCCESS);
2380: }
2382: /*@
2383: VecGetBindingPropagates - Gets whether the state of being bound to the CPU for a GPU vector type propagates to child and some other associated objects
2385: Input Parameter:
2386: . v - the vector
2388: Output Parameter:
2389: . flg - flag indicating whether the boundtocpu flag will be propagated
2391: Level: developer
2393: .seealso: [](ch_vectors), `Vec`, `VecSetBindingPropagates()`
2394: @*/
2395: PetscErrorCode VecGetBindingPropagates(Vec v, PetscBool *flg)
2396: {
2397: PetscFunctionBegin;
2399: PetscAssertPointer(flg, 2);
2400: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2401: *flg = v->bindingpropagates;
2402: #else
2403: *flg = PETSC_FALSE;
2404: #endif
2405: PetscFunctionReturn(PETSC_SUCCESS);
2406: }
2408: /*@C
2409: VecSetPinnedMemoryMin - Set the minimum data size for which pinned memory will be used for host (CPU) allocations.
2411: Logically Collective
2413: Input Parameters:
2414: + v - the vector
2415: - mbytes - minimum data size in bytes
2417: Options Database Key:
2418: . -vec_pinned_memory_min size - minimum size (in bytes) for an allocation to use pinned memory on host.
2420: Level: developer
2422: Note:
2423: Specifying -1 ensures that pinned memory will never be used.
2425: .seealso: [](ch_vectors), `Vec`, `VecGetPinnedMemoryMin()`
2426: @*/
2427: PetscErrorCode VecSetPinnedMemoryMin(Vec v, size_t mbytes)
2428: {
2429: PetscFunctionBegin;
2431: #if PetscDefined(HAVE_DEVICE)
2432: v->minimum_bytes_pinned_memory = mbytes;
2433: #endif
2434: PetscFunctionReturn(PETSC_SUCCESS);
2435: }
2437: /*@C
2438: VecGetPinnedMemoryMin - Get the minimum data size for which pinned memory will be used for host (CPU) allocations.
2440: Logically Collective
2442: Input Parameter:
2443: . v - the vector
2445: Output Parameter:
2446: . mbytes - minimum data size in bytes
2448: Level: developer
2450: .seealso: [](ch_vectors), `Vec`, `VecSetPinnedMemoryMin()`
2451: @*/
2452: PetscErrorCode VecGetPinnedMemoryMin(Vec v, size_t *mbytes)
2453: {
2454: PetscFunctionBegin;
2456: PetscAssertPointer(mbytes, 2);
2457: #if PetscDefined(HAVE_DEVICE)
2458: *mbytes = v->minimum_bytes_pinned_memory;
2459: #endif
2460: PetscFunctionReturn(PETSC_SUCCESS);
2461: }
2463: /*@
2464: VecGetOffloadMask - Get the offload mask of a `Vec`
2466: Not Collective
2468: Input Parameter:
2469: . v - the vector
2471: Output Parameter:
2472: . mask - corresponding `PetscOffloadMask` enum value.
2474: Level: intermediate
2476: .seealso: [](ch_vectors), `Vec`, `VecCreateSeqCUDA()`, `VecCreateSeqViennaCL()`, `VecGetArray()`, `VecGetType()`
2477: @*/
2478: PetscErrorCode VecGetOffloadMask(Vec v, PetscOffloadMask *mask)
2479: {
2480: PetscFunctionBegin;
2482: PetscAssertPointer(mask, 2);
2483: *mask = v->offloadmask;
2484: PetscFunctionReturn(PETSC_SUCCESS);
2485: }
2487: #if !defined(PETSC_HAVE_VIENNACL)
2488: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLContext(Vec v, PETSC_UINTPTR_T *ctx)
2489: {
2490: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_context");
2491: }
2493: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLQueue(Vec v, PETSC_UINTPTR_T *queue)
2494: {
2495: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_command_queue");
2496: }
2498: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMem(Vec v, PETSC_UINTPTR_T *queue)
2499: {
2500: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2501: }
2503: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMemRead(Vec v, PETSC_UINTPTR_T *queue)
2504: {
2505: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2506: }
2508: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMemWrite(Vec v, PETSC_UINTPTR_T *queue)
2509: {
2510: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2511: }
2513: PETSC_EXTERN PetscErrorCode VecViennaCLRestoreCLMemWrite(Vec v)
2514: {
2515: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to restore a Vec's cl_mem");
2516: }
2517: #endif
2519: 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)
2520: {
2521: const PetscScalar *u, *y;
2522: const PetscScalar *atola = NULL, *rtola = NULL, *erra = NULL;
2523: PetscInt n, n_loc = 0, na_loc = 0, nr_loc = 0;
2524: PetscReal nrm = 0, nrma = 0, nrmr = 0, err_loc[6];
2526: PetscFunctionBegin;
2527: #define SkipSmallValue(a, b, tol) \
2528: if (PetscAbsScalar(a) < tol || PetscAbsScalar(b) < tol) continue
2530: PetscCall(VecGetLocalSize(U, &n));
2531: PetscCall(VecGetArrayRead(U, &u));
2532: PetscCall(VecGetArrayRead(Y, &y));
2533: if (E) PetscCall(VecGetArrayRead(E, &erra));
2534: if (vatol) PetscCall(VecGetArrayRead(vatol, &atola));
2535: if (vrtol) PetscCall(VecGetArrayRead(vrtol, &rtola));
2536: for (PetscInt i = 0; i < n; i++) {
2537: PetscReal err, tol, tola, tolr;
2539: SkipSmallValue(y[i], u[i], ignore_max);
2540: atol = atola ? PetscRealPart(atola[i]) : atol;
2541: rtol = rtola ? PetscRealPart(rtola[i]) : rtol;
2542: err = erra ? PetscAbsScalar(erra[i]) : PetscAbsScalar(y[i] - u[i]);
2543: tola = atol;
2544: tolr = rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
2545: tol = tola + tolr;
2546: if (tola > 0.) {
2547: if (wnormtype == NORM_INFINITY) nrma = PetscMax(nrma, err / tola);
2548: else nrma += PetscSqr(err / tola);
2549: na_loc++;
2550: }
2551: if (tolr > 0.) {
2552: if (wnormtype == NORM_INFINITY) nrmr = PetscMax(nrmr, err / tolr);
2553: else nrmr += PetscSqr(err / tolr);
2554: nr_loc++;
2555: }
2556: if (tol > 0.) {
2557: if (wnormtype == NORM_INFINITY) nrm = PetscMax(nrm, err / tol);
2558: else nrm += PetscSqr(err / tol);
2559: n_loc++;
2560: }
2561: }
2562: if (E) PetscCall(VecRestoreArrayRead(E, &erra));
2563: if (vatol) PetscCall(VecRestoreArrayRead(vatol, &atola));
2564: if (vrtol) PetscCall(VecRestoreArrayRead(vrtol, &rtola));
2565: PetscCall(VecRestoreArrayRead(U, &u));
2566: PetscCall(VecRestoreArrayRead(Y, &y));
2567: #undef SkipSmallValue
2569: err_loc[0] = nrm;
2570: err_loc[1] = nrma;
2571: err_loc[2] = nrmr;
2572: err_loc[3] = (PetscReal)n_loc;
2573: err_loc[4] = (PetscReal)na_loc;
2574: err_loc[5] = (PetscReal)nr_loc;
2575: if (wnormtype == NORM_2) {
2576: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, err_loc, 6, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)U)));
2577: } else {
2578: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, err_loc, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)U)));
2579: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, err_loc + 3, 3, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)U)));
2580: }
2581: if (wnormtype == NORM_2) {
2582: *norm = PetscSqrtReal(err_loc[0]);
2583: *norma = PetscSqrtReal(err_loc[1]);
2584: *normr = PetscSqrtReal(err_loc[2]);
2585: } else {
2586: *norm = err_loc[0];
2587: *norma = err_loc[1];
2588: *normr = err_loc[2];
2589: }
2590: *norm_loc = (PetscInt)err_loc[3];
2591: *norma_loc = (PetscInt)err_loc[4];
2592: *normr_loc = (PetscInt)err_loc[5];
2593: PetscFunctionReturn(PETSC_SUCCESS);
2594: }
2596: /*@
2597: VecErrorWeightedNorms - compute a weighted norm of the difference between two vectors
2599: Collective
2601: Input Parameters:
2602: + U - first vector to be compared
2603: . Y - second vector to be compared
2604: . E - optional third vector representing the error (if not provided, the error is ||U-Y||)
2605: . wnormtype - norm type
2606: . atol - scalar for absolute tolerance
2607: . vatol - vector representing per-entry absolute tolerances (can be ``NULL``)
2608: . rtol - scalar for relative tolerance
2609: . vrtol - vector representing per-entry relative tolerances (can be ``NULL``)
2610: - ignore_max - ignore values smaller than this value in absolute terms.
2612: Output Parameters:
2613: + norm - weighted norm
2614: . norm_loc - number of vector locations used for the weighted norm
2615: . norma - weighted norm based on the absolute tolerance
2616: . norma_loc - number of vector locations used for the absolute weighted norm
2617: . normr - weighted norm based on the relative tolerance
2618: - normr_loc - number of vector locations used for the relative weighted norm
2620: Level: developer
2622: Notes:
2623: This is primarily used for computing weighted local truncation errors in ``TS``.
2625: .seealso: [](ch_vectors), `Vec`, `NormType`, `TSErrorWeightedNorm()`, `TSErrorWeightedENorm()`
2626: @*/
2627: 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)
2628: {
2629: PetscFunctionBegin;
2634: if (E) {
2637: }
2640: if (vatol) {
2643: }
2645: if (vrtol) {
2648: }
2650: PetscAssertPointer(norm, 10);
2651: PetscAssertPointer(norm_loc, 11);
2652: PetscAssertPointer(norma, 12);
2653: PetscAssertPointer(norma_loc, 13);
2654: PetscAssertPointer(normr, 14);
2655: PetscAssertPointer(normr_loc, 15);
2656: PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)U), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
2658: /* There are potentially 5 vectors involved, some of them may happen to be of different type or bound to cpu.
2659: Here we check that they all implement the same operation and call it if so.
2660: Otherwise, we call the _Basic implementation that always works (provided VecGetArrayRead is implemented). */
2661: PetscBool sameop = (PetscBool)(U->ops->errorwnorm && U->ops->errorwnorm == Y->ops->errorwnorm);
2662: if (sameop && E) sameop = (PetscBool)(U->ops->errorwnorm == E->ops->errorwnorm);
2663: if (sameop && vatol) sameop = (PetscBool)(U->ops->errorwnorm == vatol->ops->errorwnorm);
2664: if (sameop && vrtol) sameop = (PetscBool)(U->ops->errorwnorm == vrtol->ops->errorwnorm);
2665: if (sameop) PetscUseTypeMethod(U, errorwnorm, Y, E, wnormtype, atol, vatol, rtol, vrtol, ignore_max, norm, norm_loc, norma, norma_loc, normr, normr_loc);
2666: else PetscCall(VecErrorWeightedNorms_Basic(U, Y, E, wnormtype, atol, vatol, rtol, vrtol, ignore_max, norm, norm_loc, norma, norma_loc, normr, normr_loc));
2667: PetscFunctionReturn(PETSC_SUCCESS);
2668: }