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_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: IS is_coo_i;
200: /* The default implementation only supports ncoo within limit of PetscInt */
201: PetscCheck(ncoo <= PETSC_MAX_INT, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncoo %" PetscCount_FMT " overflowed PetscInt; configure --with-64-bit-indices or request support", ncoo);
202: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncoo, coo_i, PETSC_COPY_VALUES, &is_coo_i));
203: PetscCall(PetscObjectCompose((PetscObject)x, "__PETSc_coo_i", (PetscObject)is_coo_i));
204: PetscCall(ISDestroy(&is_coo_i));
205: }
206: PetscCall(PetscLogEventEnd(VEC_SetPreallocateCOO, x, 0, 0, 0));
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: /*@
211: VecSetPreallocationCOOLocal - set preallocation for vectors using a coordinate format of the entries with local indices
213: Collective
215: Input Parameters:
216: + x - vector being preallocated
217: . ncoo - number of entries
218: - coo_i - row indices (local numbering; may be modified)
220: Level: beginner
222: Notes:
223: This and `VecSetValuesCOO()` provide an alternative API to using `VecSetValuesLocal()` to provide vector values.
225: This API is particularly efficient for use on GPUs.
227: The local indices are translated using the local to global mapping, thus `VecSetLocalToGlobalMapping()` must have been
228: called prior to this function.
230: The indices coo_i may be modified within this function. They might be translated to corresponding global
231: indices, but the caller should not rely on them having any specific value after this function returns. The arrays
232: can be freed or reused immediately after this function returns.
234: Entries can be repeated. Negative indices and remote indices might be allowed. see `VecSetPreallocationCOO()`.
236: .seealso: [](ch_vectors), `Vec`, `VecSetPreallocationCOO()`, `VecSetValuesCOO()`
237: @*/
238: PetscErrorCode VecSetPreallocationCOOLocal(Vec x, PetscCount ncoo, PetscInt coo_i[])
239: {
240: ISLocalToGlobalMapping ltog;
242: PetscFunctionBegin;
245: if (ncoo) PetscAssertPointer(coo_i, 3);
246: PetscCheck(ncoo <= PETSC_MAX_INT, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncoo %" PetscCount_FMT " overflowed PetscInt; configure --with-64-bit-indices or request support", ncoo);
247: PetscCall(PetscLayoutSetUp(x->map));
248: PetscCall(VecGetLocalToGlobalMapping(x, <og));
249: if (ltog) PetscCall(ISLocalToGlobalMappingApply(ltog, ncoo, coo_i, coo_i));
250: PetscCall(VecSetPreallocationCOO(x, ncoo, coo_i));
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: /*@
255: VecSetValuesCOO - set values at once in a vector preallocated using `VecSetPreallocationCOO()`
257: Collective
259: Input Parameters:
260: + x - vector being set
261: . coo_v - the value array
262: - imode - the insert mode
264: Level: beginner
266: Note:
267: This and `VecSetPreallocationCOO() or ``VecSetPreallocationCOOLocal()` provide an alternative API to using `VecSetValues()` to provide vector values.
269: This API is particularly efficient for use on GPUs.
271: The values must follow the order of the indices prescribed with `VecSetPreallocationCOO()` or `VecSetPreallocationCOOLocal()`.
272: When repeated entries are specified in the COO indices the `coo_v` values are first properly summed, regardless of the value of `imode`.
273: The imode flag indicates if `coo_v` must be added to the current values of the vector (`ADD_VALUES`) or overwritten (`INSERT_VALUES`).
274: `VecAssemblyBegin()` and `VecAssemblyEnd()` do not need to be called after this routine. It automatically handles the assembly process.
276: .seealso: [](ch_vectors), `Vec`, `VecSetPreallocationCOO()`, `VecSetPreallocationCOOLocal()`, `VecSetValues()`
277: @*/
278: PetscErrorCode VecSetValuesCOO(Vec x, const PetscScalar coo_v[], InsertMode imode)
279: {
280: PetscFunctionBegin;
284: PetscCall(PetscLogEventBegin(VEC_SetValuesCOO, x, 0, 0, 0));
285: if (x->ops->setvaluescoo) {
286: PetscUseTypeMethod(x, setvaluescoo, coo_v, imode);
287: PetscCall(PetscObjectStateIncrease((PetscObject)x));
288: } else {
289: IS is_coo_i;
290: const PetscInt *coo_i;
291: PetscInt ncoo;
292: PetscMemType mtype;
294: PetscCall(PetscGetMemType(coo_v, &mtype));
295: PetscCheck(mtype == PETSC_MEMTYPE_HOST, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_WRONG, "The basic VecSetValuesCOO() only supports v[] on host");
296: PetscCall(PetscObjectQuery((PetscObject)x, "__PETSc_coo_i", (PetscObject *)&is_coo_i));
297: PetscCheck(is_coo_i, PetscObjectComm((PetscObject)x), PETSC_ERR_COR, "Missing coo_i IS");
298: PetscCall(ISGetLocalSize(is_coo_i, &ncoo));
299: PetscCall(ISGetIndices(is_coo_i, &coo_i));
300: if (imode != ADD_VALUES) PetscCall(VecZeroEntries(x));
301: PetscCall(VecSetValues(x, ncoo, coo_i, coo_v, ADD_VALUES));
302: PetscCall(ISRestoreIndices(is_coo_i, &coo_i));
303: PetscCall(VecAssemblyBegin(x));
304: PetscCall(VecAssemblyEnd(x));
305: }
306: PetscCall(PetscLogEventEnd(VEC_SetValuesCOO, x, 0, 0, 0));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: 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))
311: {
312: PetscErrorCode (*async_fn)(Vec, Vec, Vec, PetscDeviceContext) = NULL;
314: PetscFunctionBegin;
321: PetscCheckSameTypeAndComm(x, 2, y, 3);
322: PetscCheckSameTypeAndComm(y, 3, w, 1);
323: VecCheckSameSize(w, 1, x, 2);
324: VecCheckSameSize(w, 1, y, 3);
325: VecCheckAssembled(x);
326: VecCheckAssembled(y);
327: PetscCall(VecSetErrorIfLocked(w, 1));
330: if (dctx) PetscCall(PetscObjectQueryFunction((PetscObject)w, async_name, &async_fn));
331: if (event) PetscCall(PetscLogEventBegin(event, x, y, w, 0));
332: if (async_fn) {
333: PetscCall((*async_fn)(w, x, y, dctx));
334: } else {
335: PetscCall((*pointwise_op)(w, x, y));
336: }
337: if (event) PetscCall(PetscLogEventEnd(event, x, y, w, 0));
338: PetscCall(PetscObjectStateIncrease((PetscObject)w));
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: PetscErrorCode VecPointwiseMaxAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
343: {
344: PetscFunctionBegin;
345: // REVIEW ME: no log event?
346: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMax), w->ops->pointwisemax));
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: /*@
351: VecPointwiseMax - Computes the component-wise maximum `w[i] = max(x[i], y[i])`.
353: Logically Collective
355: Input Parameters:
356: + x - the first input vector
357: - y - the second input vector
359: Output Parameter:
360: . w - the result
362: Level: advanced
364: Notes:
365: Any subset of the `x`, `y`, and `w` may be the same vector.
367: For complex numbers compares only the real part
369: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
370: @*/
371: PetscErrorCode VecPointwiseMax(Vec w, Vec x, Vec y)
372: {
373: PetscFunctionBegin;
374: PetscCall(VecPointwiseMaxAsync_Private(w, x, y, NULL));
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: PetscErrorCode VecPointwiseMinAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
379: {
380: PetscFunctionBegin;
381: // REVIEW ME: no log event?
382: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMin), w->ops->pointwisemin));
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: /*@
387: VecPointwiseMin - Computes the component-wise minimum `w[i] = min(x[i], y[i])`.
389: Logically Collective
391: Input Parameters:
392: + x - the first input vector
393: - y - the second input vector
395: Output Parameter:
396: . w - the result
398: Level: advanced
400: Notes:
401: Any subset of the `x`, `y`, and `w` may be the same vector.
403: For complex numbers compares only the real part
405: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
406: @*/
407: PetscErrorCode VecPointwiseMin(Vec w, Vec x, Vec y)
408: {
409: PetscFunctionBegin;
410: PetscCall(VecPointwiseMinAsync_Private(w, x, y, NULL));
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: PetscErrorCode VecPointwiseMaxAbsAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
415: {
416: PetscFunctionBegin;
417: // REVIEW ME: no log event?
418: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMaxAbs), w->ops->pointwisemaxabs));
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: /*@
423: VecPointwiseMaxAbs - Computes the component-wise maximum of the absolute values `w[i] = max(abs(x[i]), abs(y[i]))`.
425: Logically Collective
427: Input Parameters:
428: + x - the first input vector
429: - y - the second input vector
431: Output Parameter:
432: . w - the result
434: Level: advanced
436: Notes:
437: Any subset of the `x`, `y`, and `w` may be the same vector.
439: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMin()`, `VecPointwiseMax()`, `VecMaxPointwiseDivide()`
440: @*/
441: PetscErrorCode VecPointwiseMaxAbs(Vec w, Vec x, Vec y)
442: {
443: PetscFunctionBegin;
444: PetscCall(VecPointwiseMaxAbsAsync_Private(w, x, y, NULL));
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: PetscErrorCode VecPointwiseDivideAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
449: {
450: PetscFunctionBegin;
451: // REVIEW ME: no log event?
452: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseDivide), w->ops->pointwisedivide));
453: PetscFunctionReturn(PETSC_SUCCESS);
454: }
456: /*@
457: VecPointwiseDivide - Computes the component-wise division `w[i] = x[i] / y[i]`.
459: Logically Collective
461: Input Parameters:
462: + x - the numerator vector
463: - y - the denominator vector
465: Output Parameter:
466: . w - the result
468: Level: advanced
470: Note:
471: Any subset of the `x`, `y`, and `w` may be the same vector.
473: .seealso: [](ch_vectors), `Vec`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
474: @*/
475: PetscErrorCode VecPointwiseDivide(Vec w, Vec x, Vec y)
476: {
477: PetscFunctionBegin;
478: PetscCall(VecPointwiseDivideAsync_Private(w, x, y, NULL));
479: PetscFunctionReturn(PETSC_SUCCESS);
480: }
482: PetscErrorCode VecPointwiseMultAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
483: {
484: PetscFunctionBegin;
486: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, VEC_PointwiseMult, VecAsyncFnName(PointwiseMult), w->ops->pointwisemult));
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: /*@
491: VecPointwiseMult - Computes the component-wise multiplication `w[i] = x[i] * y[i]`.
493: Logically Collective
495: Input Parameters:
496: + x - the first vector
497: - y - the second vector
499: Output Parameter:
500: . w - the result
502: Level: advanced
504: Note:
505: Any subset of the `x`, `y`, and `w` may be the same vector.
507: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
508: @*/
509: PetscErrorCode VecPointwiseMult(Vec w, Vec x, Vec y)
510: {
511: PetscFunctionBegin;
512: PetscCall(VecPointwiseMultAsync_Private(w, x, y, NULL));
513: PetscFunctionReturn(PETSC_SUCCESS);
514: }
516: /*@
517: VecDuplicate - Creates a new vector of the same type as an existing vector.
519: Collective
521: Input Parameter:
522: . v - a vector to mimic
524: Output Parameter:
525: . newv - location to put new vector
527: Level: beginner
529: Notes:
530: `VecDuplicate()` DOES NOT COPY the vector entries, but rather allocates storage
531: for the new vector. Use `VecCopy()` to copy a vector.
533: Use `VecDestroy()` to free the space. Use `VecDuplicateVecs()` to get several
534: vectors.
536: .seealso: [](ch_vectors), `Vec`, `VecDestroy()`, `VecDuplicateVecs()`, `VecCreate()`, `VecCopy()`
537: @*/
538: PetscErrorCode VecDuplicate(Vec v, Vec *newv)
539: {
540: PetscFunctionBegin;
542: PetscAssertPointer(newv, 2);
544: PetscUseTypeMethod(v, duplicate, newv);
545: #if PetscDefined(HAVE_DEVICE)
546: if (v->boundtocpu && v->bindingpropagates) {
547: PetscCall(VecSetBindingPropagates(*newv, PETSC_TRUE));
548: PetscCall(VecBindToCPU(*newv, PETSC_TRUE));
549: }
550: #endif
551: PetscCall(PetscObjectStateIncrease((PetscObject)*newv));
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
555: /*@C
556: VecDestroy - Destroys a vector.
558: Collective
560: Input Parameter:
561: . v - the vector
563: Level: beginner
565: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecDuplicate()`, `VecDestroyVecs()`
566: @*/
567: PetscErrorCode VecDestroy(Vec *v)
568: {
569: PetscFunctionBegin;
570: PetscAssertPointer(v, 1);
571: if (!*v) PetscFunctionReturn(PETSC_SUCCESS);
573: if (--((PetscObject)*v)->refct > 0) {
574: *v = NULL;
575: PetscFunctionReturn(PETSC_SUCCESS);
576: }
578: PetscCall(PetscObjectSAWsViewOff((PetscObject)*v));
579: /* destroy the internal part */
580: PetscTryTypeMethod(*v, destroy);
581: PetscCall(PetscFree((*v)->defaultrandtype));
582: /* destroy the external/common part */
583: PetscCall(PetscLayoutDestroy(&(*v)->map));
584: PetscCall(PetscHeaderDestroy(v));
585: PetscFunctionReturn(PETSC_SUCCESS);
586: }
588: /*@C
589: VecDuplicateVecs - Creates several vectors of the same type as an existing vector.
591: Collective
593: Input Parameters:
594: + m - the number of vectors to obtain
595: - v - a vector to mimic
597: Output Parameter:
598: . V - location to put pointer to array of vectors
600: Level: intermediate
602: Note:
603: Use `VecDestroyVecs()` to free the space. Use `VecDuplicate()` to form a single
604: vector.
606: Fortran Notes:
607: The Fortran interface is slightly different from that given below, it
608: requires one to pass in `V` a `Vec` array of size at least `m`.
609: See the [](ch_fortran) for details.
611: .seealso: [](ch_vectors), `Vec`, [](ch_fortran), `VecDestroyVecs()`, `VecDuplicate()`, `VecCreate()`, `VecDuplicateVecsF90()`
612: @*/
613: PetscErrorCode VecDuplicateVecs(Vec v, PetscInt m, Vec *V[])
614: {
615: PetscFunctionBegin;
617: PetscAssertPointer(V, 3);
619: PetscUseTypeMethod(v, duplicatevecs, m, V);
620: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
621: if (v->boundtocpu && v->bindingpropagates) {
622: PetscInt i;
624: for (i = 0; i < m; i++) {
625: /* Since ops->duplicatevecs might itself propagate the value of boundtocpu,
626: * avoid unnecessary overhead by only calling VecBindToCPU() if the vector isn't already bound. */
627: if (!(*V)[i]->boundtocpu) {
628: PetscCall(VecSetBindingPropagates((*V)[i], PETSC_TRUE));
629: PetscCall(VecBindToCPU((*V)[i], PETSC_TRUE));
630: }
631: }
632: }
633: #endif
634: PetscFunctionReturn(PETSC_SUCCESS);
635: }
637: /*@C
638: VecDestroyVecs - Frees a block of vectors obtained with `VecDuplicateVecs()`.
640: Collective
642: Input Parameters:
643: + m - the number of vectors previously obtained, if zero no vectors are destroyed
644: - vv - pointer to pointer to array of vector pointers, if `NULL` no vectors are destroyed
646: Level: intermediate
648: Fortran Notes:
649: The Fortran interface is slightly different from that given below.
650: See the [](ch_fortran) for details.
652: .seealso: [](ch_vectors), `Vec`, [](ch_fortran), `VecDuplicateVecs()`, `VecDestroyVecsf90()`
653: @*/
654: PetscErrorCode VecDestroyVecs(PetscInt m, Vec *vv[])
655: {
656: PetscFunctionBegin;
657: PetscAssertPointer(vv, 2);
658: PetscCheck(m >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Trying to destroy negative number of vectors %" PetscInt_FMT, m);
659: if (!m || !*vv) {
660: *vv = NULL;
661: PetscFunctionReturn(PETSC_SUCCESS);
662: }
665: PetscCall((*(**vv)->ops->destroyvecs)(m, *vv));
666: *vv = NULL;
667: PetscFunctionReturn(PETSC_SUCCESS);
668: }
670: /*@C
671: VecViewFromOptions - View a vector based on values in the options database
673: Collective
675: Input Parameters:
676: + A - the vector
677: . obj - Optional object that provides the options prefix for this viewing
678: - name - command line option
680: Level: intermediate
682: Note:
683: See `PetscObjectViewFromOptions()` to see the `PetscViewer` and PetscViewerFormat` available
685: .seealso: [](ch_vectors), `Vec`, `VecView`, `PetscObjectViewFromOptions()`, `VecCreate()`
686: @*/
687: PetscErrorCode VecViewFromOptions(Vec A, PetscObject obj, const char name[])
688: {
689: PetscFunctionBegin;
691: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: /*@C
696: VecView - Views a vector object.
698: Collective
700: Input Parameters:
701: + vec - the vector
702: - viewer - an optional `PetscViewer` visualization context
704: Level: beginner
706: Notes:
707: The available visualization contexts include
708: + `PETSC_VIEWER_STDOUT_SELF` - for sequential vectors
709: . `PETSC_VIEWER_STDOUT_WORLD` - for parallel vectors created on `PETSC_COMM_WORLD`
710: - `PETSC_VIEWER_STDOUT`_(comm) - for parallel vectors created on MPI communicator comm
712: You can change the format the vector is printed using the
713: option `PetscViewerPushFormat()`.
715: The user can open alternative viewers with
716: + `PetscViewerASCIIOpen()` - Outputs vector to a specified file
717: . `PetscViewerBinaryOpen()` - Outputs vector in binary to a
718: specified file; corresponding input uses `VecLoad()`
719: . `PetscViewerDrawOpen()` - Outputs vector to an X window display
720: . `PetscViewerSocketOpen()` - Outputs vector to Socket viewer
721: - `PetscViewerHDF5Open()` - Outputs vector to HDF5 file viewer
723: The user can call `PetscViewerPushFormat()` to specify the output
724: format of ASCII printed objects (when using `PETSC_VIEWER_STDOUT_SELF`,
725: `PETSC_VIEWER_STDOUT_WORLD` and `PetscViewerASCIIOpen()`). Available formats include
726: + `PETSC_VIEWER_DEFAULT` - default, prints vector contents
727: . `PETSC_VIEWER_ASCII_MATLAB` - prints vector contents in MATLAB format
728: . `PETSC_VIEWER_ASCII_INDEX` - prints vector contents, including indices of vector elements
729: - `PETSC_VIEWER_ASCII_COMMON` - prints vector contents, using a
730: format common among all vector types
732: You can pass any number of vector objects, or other PETSc objects to the same viewer.
734: In the debugger you can do call `VecView`(v,0) to display the vector. (The same holds for any PETSc object viewer).
736: Notes for binary viewer:
737: If you pass multiple vectors to a binary viewer you can read them back in the same order
738: with `VecLoad()`.
740: If the blocksize of the vector is greater than one then you must provide a unique prefix to
741: the vector with `PetscObjectSetOptionsPrefix`((`PetscObject`)vec,"uniqueprefix"); BEFORE calling `VecView()` on the
742: vector to be stored and then set that same unique prefix on the vector that you pass to `VecLoad()`. The blocksize
743: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
744: filename. If you copy the binary file, make sure you copy the associated .info file with it.
746: See the manual page for `VecLoad()` on the exact format the binary viewer stores
747: the values in the file.
749: Notes for HDF5 Viewer:
750: The name of the `Vec` (given with `PetscObjectSetName()` is the name that is used
751: for the object in the HDF5 file. If you wish to store the same Vec into multiple
752: datasets in the same file (typically with different values), you must change its
753: name each time before calling the `VecView()`. To load the same vector,
754: the name of the Vec object passed to `VecLoad()` must be the same.
756: If the block size of the vector is greater than 1 then it is used as the first dimension in the HDF5 array.
757: If the function `PetscViewerHDF5SetBaseDimension2()`is called then even if the block size is one it will
758: be used as the first dimension in the HDF5 array (that is the HDF5 array will always be two dimensional)
759: See also `PetscViewerHDF5SetTimestep()` which adds an additional complication to reading and writing `Vec`
760: with the HDF5 viewer.
762: .seealso: [](ch_vectors), `Vec`, `VecViewFromOptions()`, `PetscViewerASCIIOpen()`, `PetscViewerDrawOpen()`, `PetscDrawLGCreate()`,
763: `PetscViewerSocketOpen()`, `PetscViewerBinaryOpen()`, `VecLoad()`, `PetscViewerCreate()`,
764: `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`, `PetscViewerHDF5SetTimestep()`
765: @*/
766: PetscErrorCode VecView(Vec vec, PetscViewer viewer)
767: {
768: PetscBool iascii;
769: PetscViewerFormat format;
770: PetscMPIInt size;
772: PetscFunctionBegin;
775: VecCheckAssembled(vec);
776: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec), &viewer));
778: PetscCall(PetscViewerGetFormat(viewer, &format));
779: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
780: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) PetscFunctionReturn(PETSC_SUCCESS);
782: PetscCheck(!vec->stash.n && !vec->bstash.n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call VecAssemblyBegin/End() before viewing this vector");
784: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
785: if (iascii) {
786: PetscInt rows, bs;
788: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)vec, viewer));
789: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
790: PetscCall(PetscViewerASCIIPushTab(viewer));
791: PetscCall(VecGetSize(vec, &rows));
792: PetscCall(VecGetBlockSize(vec, &bs));
793: if (bs != 1) {
794: PetscCall(PetscViewerASCIIPrintf(viewer, "length=%" PetscInt_FMT ", bs=%" PetscInt_FMT "\n", rows, bs));
795: } else {
796: PetscCall(PetscViewerASCIIPrintf(viewer, "length=%" PetscInt_FMT "\n", rows));
797: }
798: PetscCall(PetscViewerASCIIPopTab(viewer));
799: }
800: }
801: PetscCall(VecLockReadPush(vec));
802: PetscCall(PetscLogEventBegin(VEC_View, vec, viewer, 0, 0));
803: if ((format == PETSC_VIEWER_NATIVE || format == PETSC_VIEWER_LOAD_BALANCE) && vec->ops->viewnative) {
804: PetscUseTypeMethod(vec, viewnative, viewer);
805: } else {
806: PetscUseTypeMethod(vec, view, viewer);
807: }
808: PetscCall(VecLockReadPop(vec));
809: PetscCall(PetscLogEventEnd(VEC_View, vec, viewer, 0, 0));
810: PetscFunctionReturn(PETSC_SUCCESS);
811: }
813: #if defined(PETSC_USE_DEBUG)
814: #include <../src/sys/totalview/tv_data_display.h>
815: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
816: {
817: const PetscScalar *values;
818: char type[32];
820: TV_add_row("Local rows", "int", &v->map->n);
821: TV_add_row("Global rows", "int", &v->map->N);
822: TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
823: PetscCall(VecGetArrayRead((Vec)v, &values));
824: PetscCall(PetscSNPrintf(type, 32, "double[%" PetscInt_FMT "]", v->map->n));
825: TV_add_row("values", type, values);
826: PetscCall(VecRestoreArrayRead((Vec)v, &values));
827: return TV_format_OK;
828: }
829: #endif
831: /*@C
832: VecViewNative - Views a vector object with the original type specific viewer
834: Collective
836: Input Parameters:
837: + vec - the vector
838: - viewer - an optional `PetscViewer` visualization context
840: Level: developer
842: Note:
843: This can be used with, for example, vectors obtained with `DMCreateGlobalVector()` for a `DMDA` to display the vector
844: in the PETSc storage format (each MPI process values follow the previous MPI processes) instead of the "natural" grid
845: ordering.
847: .seealso: [](ch_vectors), `Vec`, `PetscViewerASCIIOpen()`, `PetscViewerDrawOpen()`, `PetscDrawLGCreate()`, `VecView()`
848: `PetscViewerSocketOpen()`, `PetscViewerBinaryOpen()`, `VecLoad()`, `PetscViewerCreate()`,
849: `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`, `PetscViewerHDF5SetTimestep()`
850: @*/
851: PetscErrorCode VecViewNative(Vec vec, PetscViewer viewer)
852: {
853: PetscFunctionBegin;
856: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec), &viewer));
858: PetscUseTypeMethod(vec, viewnative, viewer);
859: PetscFunctionReturn(PETSC_SUCCESS);
860: }
862: /*@
863: VecGetSize - Returns the global number of elements of the vector.
865: Not Collective
867: Input Parameter:
868: . x - the vector
870: Output Parameter:
871: . size - the global length of the vector
873: Level: beginner
875: .seealso: [](ch_vectors), `Vec`, `VecGetLocalSize()`
876: @*/
877: PetscErrorCode VecGetSize(Vec x, PetscInt *size)
878: {
879: PetscFunctionBegin;
881: PetscAssertPointer(size, 2);
883: PetscUseTypeMethod(x, getsize, size);
884: PetscFunctionReturn(PETSC_SUCCESS);
885: }
887: /*@
888: VecGetLocalSize - Returns the number of elements of the vector stored
889: in local memory (that is on this MPI process)
891: Not Collective
893: Input Parameter:
894: . x - the vector
896: Output Parameter:
897: . size - the length of the local piece of the vector
899: Level: beginner
901: .seealso: [](ch_vectors), `Vec`, `VecGetSize()`
902: @*/
903: PetscErrorCode VecGetLocalSize(Vec x, PetscInt *size)
904: {
905: PetscFunctionBegin;
907: PetscAssertPointer(size, 2);
909: PetscUseTypeMethod(x, getlocalsize, size);
910: PetscFunctionReturn(PETSC_SUCCESS);
911: }
913: /*@C
914: VecGetOwnershipRange - Returns the range of indices owned by
915: this process. The vector is laid out with the
916: first `n1` elements on the first processor, next `n2` elements on the
917: second, etc. For certain parallel layouts this range may not be
918: well defined.
920: Not Collective
922: Input Parameter:
923: . x - the vector
925: Output Parameters:
926: + low - the first local element, pass in `NULL` if not interested
927: - high - one more than the last local element, pass in `NULL` if not interested
929: Level: beginner
931: Notes:
932: If the `Vec` was obtained from a `DM` with `DMCreateGlobalVector()`, then the range values are determined by the specific `DM`.
934: If the `Vec` was created directly the range values are determined by the local size passed to `VecSetSizes()` or `VecCreateMPI()`.
935: If `PETSC_DECIDE` was passed as the local size, then the vector uses default values for the range using `PetscSplitOwnership()`.
937: The high argument is one more than the last element stored locally.
939: For certain `DM`, such as `DMDA`, it is better to use `DM` specific routines, such as `DMDAGetGhostCorners()`, to determine
940: the local values in the vector.
942: .seealso: [](ch_vectors), `Vec`, `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `VecGetOwnershipRanges()`, `PetscSplitOwnership()`,
943: `VecSetSizes()`, `VecCreateMPI()`, `PetscLayout`, `DMDAGetGhostCorners()`, `DM`
944: @*/
945: PetscErrorCode VecGetOwnershipRange(Vec x, PetscInt *low, PetscInt *high)
946: {
947: PetscFunctionBegin;
950: if (low) PetscAssertPointer(low, 2);
951: if (high) PetscAssertPointer(high, 3);
952: if (low) *low = x->map->rstart;
953: if (high) *high = x->map->rend;
954: PetscFunctionReturn(PETSC_SUCCESS);
955: }
957: /*@C
958: VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
959: The vector is laid out with the
960: first `n1` elements on the first processor, next `n2` elements on the
961: second, etc. For certain parallel layouts this range may not be
962: well defined.
964: Not Collective
966: Input Parameter:
967: . x - the vector
969: Output Parameter:
970: . ranges - array of length `size` + 1 with the start and end+1 for each process
972: Level: beginner
974: Notes:
975: If the `Vec` was obtained from a `DM` with `DMCreateGlobalVector()`, then the range values are determined by the specific `DM`.
977: If the `Vec` was created directly the range values are determined by the local size passed to `VecSetSizes()` or `VecCreateMPI()`.
978: If `PETSC_DECIDE` was passed as the local size, then the vector uses default values for the range using `PetscSplitOwnership()`.
980: The high argument is one more than the last element stored locally.
982: For certain `DM`, such as `DMDA`, it is better to use `DM` specific routines, such as `DMDAGetGhostCorners()`, to determine
983: the local values in the vector.
985: The high argument is one more than the last element stored locally.
987: If `ranges` are used after all vectors that share the ranges has been destroyed, then the program will crash accessing `ranges`.
989: Fortran Notes:
990: You must PASS in an array of length `size` + 1, where `size` is the size of the communicator owning the vector
992: .seealso: [](ch_vectors), `Vec`, `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `VecGetOwnershipRange()`, `PetscSplitOwnership()`,
993: `VecSetSizes()`, `VecCreateMPI()`, `PetscLayout`, `DMDAGetGhostCorners()`, `DM`
994: @*/
995: PetscErrorCode VecGetOwnershipRanges(Vec x, const PetscInt *ranges[])
996: {
997: PetscFunctionBegin;
1000: PetscCall(PetscLayoutGetRanges(x->map, ranges));
1001: PetscFunctionReturn(PETSC_SUCCESS);
1002: }
1004: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
1005: /*@
1006: VecSetOption - Sets an option for controlling a vector's behavior.
1008: Collective
1010: Input Parameters:
1011: + x - the vector
1012: . op - the option
1013: - flag - turn the option on or off
1015: Supported Options:
1016: + `VEC_IGNORE_OFF_PROC_ENTRIES` - which causes `VecSetValues()` to ignore
1017: entries destined to be stored on a separate processor. This can be used
1018: to eliminate the global reduction in the `VecAssemblyBegin()` if you know
1019: that you have only used `VecSetValues()` to set local elements
1020: . `VEC_IGNORE_NEGATIVE_INDICES` - which means you can pass negative indices
1021: in ix in calls to `VecSetValues()` or `VecGetValues()`. These rows are simply
1022: ignored.
1023: - `VEC_SUBSET_OFF_PROC_ENTRIES` - which causes `VecAssemblyBegin()` to assume that the off-process
1024: entries will always be a subset (possibly equal) of the off-process entries set on the
1025: first assembly which had a true `VEC_SUBSET_OFF_PROC_ENTRIES` and the vector has not
1026: changed this flag afterwards. If this assembly is not such first assembly, then this
1027: assembly can reuse the communication pattern setup in that first assembly, thus avoiding
1028: a global reduction. Subsequent assemblies setting off-process values should use the same
1029: InsertMode as the first assembly.
1031: Level: intermediate
1033: Developer Notes:
1034: The `InsertMode` restriction could be removed by packing the stash messages out of place.
1036: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`
1037: @*/
1038: PetscErrorCode VecSetOption(Vec x, VecOption op, PetscBool flag)
1039: {
1040: PetscFunctionBegin;
1043: PetscTryTypeMethod(x, setoption, op, flag);
1044: PetscFunctionReturn(PETSC_SUCCESS);
1045: }
1047: /* Default routines for obtaining and releasing; */
1048: /* may be used by any implementation */
1049: PetscErrorCode VecDuplicateVecs_Default(Vec w, PetscInt m, Vec *V[])
1050: {
1051: PetscFunctionBegin;
1052: PetscCheck(m > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "m must be > 0: m = %" PetscInt_FMT, m);
1053: PetscCall(PetscMalloc1(m, V));
1054: for (PetscInt i = 0; i < m; i++) PetscCall(VecDuplicate(w, *V + i));
1055: PetscFunctionReturn(PETSC_SUCCESS);
1056: }
1058: PetscErrorCode VecDestroyVecs_Default(PetscInt m, Vec v[])
1059: {
1060: PetscInt i;
1062: PetscFunctionBegin;
1063: PetscAssertPointer(v, 2);
1064: for (i = 0; i < m; i++) PetscCall(VecDestroy(&v[i]));
1065: PetscCall(PetscFree(v));
1066: PetscFunctionReturn(PETSC_SUCCESS);
1067: }
1069: /*@
1070: VecResetArray - Resets a vector to use its default memory. Call this
1071: after the use of `VecPlaceArray()`.
1073: Not Collective
1075: Input Parameter:
1076: . vec - the vector
1078: Level: developer
1080: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecPlaceArray()`
1081: @*/
1082: PetscErrorCode VecResetArray(Vec vec)
1083: {
1084: PetscFunctionBegin;
1087: PetscUseTypeMethod(vec, resetarray);
1088: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
1089: PetscFunctionReturn(PETSC_SUCCESS);
1090: }
1092: /*@C
1093: VecLoad - Loads a vector that has been stored in binary or HDF5 format
1094: with `VecView()`.
1096: Collective
1098: Input Parameters:
1099: + vec - the newly loaded vector, this needs to have been created with `VecCreate()` or
1100: some related function before the call to `VecLoad()`.
1101: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()` or
1102: HDF5 file viewer, obtained from `PetscViewerHDF5Open()`
1104: Level: intermediate
1106: Notes:
1107: Defaults to the standard `VECSEQ` or `VECMPI`, if you want some other type of `Vec` call `VecSetFromOptions()`
1108: before calling this.
1110: The input file must contain the full global vector, as
1111: written by the routine `VecView()`.
1113: If the type or size of `vec` is not set before a call to `VecLoad()`, PETSc
1114: sets the type and the local and global sizes based on the vector it is reading in. If type and/or
1115: sizes are already set, then the same are used.
1117: If using the binary viewer and the blocksize of the vector is greater than one then you must provide a unique prefix to
1118: the vector with `PetscObjectSetOptionsPrefix`((`PetscObject`)vec,"uniqueprefix"); BEFORE calling `VecView()` on the
1119: vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
1120: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
1121: filename. If you copy the binary file, make sure you copy the associated .info file with it.
1123: If using HDF5, you must assign the `Vec` the same name as was used in the Vec
1124: that was stored in the file using `PetscObjectSetName(). Otherwise you will
1125: get the error message: "Cannot H5DOpen2() with `Vec` name NAMEOFOBJECT".
1127: If the HDF5 file contains a two dimensional array the first dimension is treated as the block size
1128: in loading the vector. Hence, for example, using MATLAB notation h5create('vector.dat','/Test_Vec',[27 1]);
1129: will load a vector of size 27 and block size 27 thus resulting in all 27 entries being on the first process of
1130: vectors communicator and the rest of the processes having zero entries
1132: Notes for advanced users when using the binary viewer:
1133: Most users should not need to know the details of the binary storage
1134: format, since `VecLoad()` and `VecView()` completely hide these details.
1135: But for anyone who's interested, the standard binary vector storage
1136: format is
1137: .vb
1138: PetscInt VEC_FILE_CLASSID
1139: PetscInt number of rows
1140: PetscScalar *values of all entries
1141: .ve
1143: In addition, PETSc automatically uses byte swapping to work on all machines; the files
1144: are written ALWAYS using big-endian ordering. On small-endian machines the numbers
1145: are converted to the small-endian format when they are read in from the file.
1146: See PetscBinaryRead() and PetscBinaryWrite() to see how this may be done.
1148: .seealso: [](ch_vectors), `Vec`, `PetscViewerBinaryOpen()`, `VecView()`, `MatLoad()`
1149: @*/
1150: PetscErrorCode VecLoad(Vec vec, PetscViewer viewer)
1151: {
1152: PetscBool isbinary, ishdf5, isadios, isexodusii;
1153: PetscViewerFormat format;
1155: PetscFunctionBegin;
1158: PetscCheckSameComm(vec, 1, viewer, 2);
1159: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1160: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1161: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERADIOS, &isadios));
1162: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWEREXODUSII, &isexodusii));
1163: PetscCheck(isbinary || ishdf5 || isadios || isexodusii, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1165: PetscCall(VecSetErrorIfLocked(vec, 1));
1166: if (!((PetscObject)vec)->type_name && !vec->ops->create) PetscCall(VecSetType(vec, VECSTANDARD));
1167: PetscCall(PetscLogEventBegin(VEC_Load, viewer, 0, 0, 0));
1168: PetscCall(PetscViewerGetFormat(viewer, &format));
1169: if (format == PETSC_VIEWER_NATIVE && vec->ops->loadnative) {
1170: PetscUseTypeMethod(vec, loadnative, viewer);
1171: } else {
1172: PetscUseTypeMethod(vec, load, viewer);
1173: }
1174: PetscCall(PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0));
1175: PetscFunctionReturn(PETSC_SUCCESS);
1176: }
1178: /*@
1179: VecReciprocal - Replaces each component of a vector by its reciprocal.
1181: Logically Collective
1183: Input Parameter:
1184: . vec - the vector
1186: Output Parameter:
1187: . vec - the vector reciprocal
1189: Level: intermediate
1191: Note:
1192: Vector entries with value 0.0 are not changed
1194: .seealso: [](ch_vectors), `Vec`, `VecLog()`, `VecExp()`, `VecSqrtAbs()`
1195: @*/
1196: PetscErrorCode VecReciprocal(Vec vec)
1197: {
1198: PetscFunctionBegin;
1199: PetscCall(VecReciprocalAsync_Private(vec, NULL));
1200: PetscFunctionReturn(PETSC_SUCCESS);
1201: }
1203: /*@C
1204: VecSetOperation - Allows the user to override a particular vector operation.
1206: Logically Collective; No Fortran Support
1208: Input Parameters:
1209: + vec - The vector to modify
1210: . op - The name of the operation
1211: - f - The function that provides the operation.
1213: Notes:
1214: `f` may be `NULL` to remove the operation from `vec`. Depending on the operation this may be
1215: allowed, however some always expect a valid function. In these cases an error will be raised
1216: when calling the interface routine in question.
1218: See `VecOperation` for an up-to-date list of override-able operations. The operations listed
1219: there have the form `VECOP_<OPERATION>`, where `<OPERATION>` is the suffix (in all capital
1220: letters) of the public interface routine (e.g., `VecView()` -> `VECOP_VIEW`).
1222: Overriding a particular `Vec`'s operation has no affect on any other `Vec`s past, present,
1223: or future. The user should also note that overriding a method is "destructive"; the previous
1224: method is not retained in any way.
1226: Level: advanced
1228: Example Usage:
1229: .vb
1230: // some new VecView() implementation, must have the same signature as the function it seeks
1231: // to replace
1232: PetscErrorCode UserVecView(Vec x, PetscViewer viewer)
1233: {
1234: PetscFunctionBeginUser;
1235: // ...
1236: PetscFunctionReturn(PETSC_SUCCESS);
1237: }
1239: // Create a VECMPI which has a pre-defined VecView() implementation
1240: VecCreateMPI(comm, n, N, &x);
1241: // Calls the VECMPI implementation for VecView()
1242: VecView(x, viewer);
1244: VecSetOperation(x, VECOP_VIEW, (void (*)(void))UserVecView);
1245: // Now calls UserVecView()
1246: VecView(x, viewer);
1247: .ve
1249: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `MatShellSetOperation()`
1250: @*/
1251: PetscErrorCode VecSetOperation(Vec vec, VecOperation op, void (*f)(void))
1252: {
1253: PetscFunctionBegin;
1255: if (op == VECOP_VIEW && !vec->ops->viewnative) {
1256: vec->ops->viewnative = vec->ops->view;
1257: } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
1258: vec->ops->loadnative = vec->ops->load;
1259: }
1260: ((void (**)(void))vec->ops)[(int)op] = f;
1261: PetscFunctionReturn(PETSC_SUCCESS);
1262: }
1264: /*@
1265: VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1266: used during the assembly process to store values that belong to
1267: other processors.
1269: Not Collective, different processes can have different size stashes
1271: Input Parameters:
1272: + vec - the vector
1273: . size - the initial size of the stash.
1274: - bsize - the initial size of the block-stash(if used).
1276: Options Database Keys:
1277: + -vecstash_initial_size <size> or <size0,size1,...sizep-1> - set initial size
1278: - -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1> - set initial block size
1280: Level: intermediate
1282: Notes:
1283: The block-stash is used for values set with `VecSetValuesBlocked()` while
1284: the stash is used for values set with `VecSetValues()`
1286: Run with the option -info and look for output of the form
1287: VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1288: to determine the appropriate value, MM, to use for size and
1289: VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1290: to determine the value, BMM to use for bsize
1292: PETSc attempts to smartly manage the stash size so there is little likelihood setting a
1293: a specific value here will affect performance
1295: .seealso: [](ch_vectors), `Vec`, `VecSetBlockSize()`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecStashView()`
1296: @*/
1297: PetscErrorCode VecStashSetInitialSize(Vec vec, PetscInt size, PetscInt bsize)
1298: {
1299: PetscFunctionBegin;
1301: PetscCall(VecStashSetInitialSize_Private(&vec->stash, size));
1302: PetscCall(VecStashSetInitialSize_Private(&vec->bstash, bsize));
1303: PetscFunctionReturn(PETSC_SUCCESS);
1304: }
1306: /*@
1307: VecSetRandom - Sets all components of a vector to random numbers.
1309: Logically Collective
1311: Input Parameters:
1312: + x - the vector
1313: - rctx - the random number context, formed by `PetscRandomCreate()`, or use `NULL` and it will create one internally.
1315: Output Parameter:
1316: . x - the vector
1318: Example of Usage:
1319: .vb
1320: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1321: VecSetRandom(x,rctx);
1322: PetscRandomDestroy(&rctx);
1323: .ve
1325: Level: intermediate
1327: .seealso: [](ch_vectors), `Vec`, `VecSet()`, `VecSetValues()`, `PetscRandomCreate()`, `PetscRandomDestroy()`
1328: @*/
1329: PetscErrorCode VecSetRandom(Vec x, PetscRandom rctx)
1330: {
1331: PetscRandom randObj = NULL;
1333: PetscFunctionBegin;
1337: VecCheckAssembled(x);
1338: PetscCall(VecSetErrorIfLocked(x, 1));
1340: if (!rctx) {
1341: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)x), &randObj));
1342: PetscCall(PetscRandomSetType(randObj, x->defaultrandtype));
1343: PetscCall(PetscRandomSetFromOptions(randObj));
1344: rctx = randObj;
1345: }
1347: PetscCall(PetscLogEventBegin(VEC_SetRandom, x, rctx, 0, 0));
1348: PetscUseTypeMethod(x, setrandom, rctx);
1349: PetscCall(PetscLogEventEnd(VEC_SetRandom, x, rctx, 0, 0));
1351: PetscCall(PetscRandomDestroy(&randObj));
1352: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1353: PetscFunctionReturn(PETSC_SUCCESS);
1354: }
1356: /*@
1357: VecZeroEntries - puts a `0.0` in each element of a vector
1359: Logically Collective
1361: Input Parameter:
1362: . vec - The vector
1364: Level: beginner
1366: Note:
1367: If the norm of the vector is known to be zero then this skips the unneeded zeroing process
1369: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecSetOptionsPrefix()`, `VecSet()`, `VecSetValues()`
1370: @*/
1371: PetscErrorCode VecZeroEntries(Vec vec)
1372: {
1373: PetscFunctionBegin;
1374: PetscCall(VecSet(vec, 0));
1375: PetscFunctionReturn(PETSC_SUCCESS);
1376: }
1378: /*
1379: VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
1380: processor and a PETSc MPI vector on more than one processor.
1382: Collective
1384: Input Parameter:
1385: . vec - The vector
1387: Level: intermediate
1389: .seealso: [](ch_vectors), `Vec`, `VecSetFromOptions()`, `VecSetType()`
1390: */
1391: static PetscErrorCode VecSetTypeFromOptions_Private(Vec vec, PetscOptionItems *PetscOptionsObject)
1392: {
1393: PetscBool opt;
1394: VecType defaultType;
1395: char typeName[256];
1396: PetscMPIInt size;
1398: PetscFunctionBegin;
1399: if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1400: else {
1401: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
1402: if (size > 1) defaultType = VECMPI;
1403: else defaultType = VECSEQ;
1404: }
1406: PetscCall(VecRegisterAll());
1407: PetscCall(PetscOptionsFList("-vec_type", "Vector type", "VecSetType", VecList, defaultType, typeName, 256, &opt));
1408: if (opt) {
1409: PetscCall(VecSetType(vec, typeName));
1410: } else {
1411: PetscCall(VecSetType(vec, defaultType));
1412: }
1413: PetscFunctionReturn(PETSC_SUCCESS);
1414: }
1416: /*@
1417: VecSetFromOptions - Configures the vector from the options database.
1419: Collective
1421: Input Parameter:
1422: . vec - The vector
1424: Level: beginner
1426: Notes:
1427: To see all options, run your program with the -help option.
1429: Must be called after `VecCreate()` but before the vector is used.
1431: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecSetOptionsPrefix()`
1432: @*/
1433: PetscErrorCode VecSetFromOptions(Vec vec)
1434: {
1435: PetscBool flg;
1436: PetscInt bind_below = 0;
1438: PetscFunctionBegin;
1441: PetscObjectOptionsBegin((PetscObject)vec);
1442: /* Handle vector type options */
1443: PetscCall(VecSetTypeFromOptions_Private(vec, PetscOptionsObject));
1445: /* Handle specific vector options */
1446: PetscTryTypeMethod(vec, setfromoptions, PetscOptionsObject);
1448: /* Bind to CPU if below a user-specified size threshold.
1449: * This perhaps belongs in the options for the GPU Vec types, but VecBindToCPU() does nothing when called on non-GPU types,
1450: * and putting it here makes is more maintainable than duplicating this for all. */
1451: 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));
1452: if (flg && vec->map->n < bind_below) PetscCall(VecBindToCPU(vec, PETSC_TRUE));
1454: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1455: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)vec, PetscOptionsObject));
1456: PetscOptionsEnd();
1457: PetscFunctionReturn(PETSC_SUCCESS);
1458: }
1460: /*@
1461: VecSetSizes - Sets the local and global sizes, and checks to determine compatibility of the sizes
1463: Collective
1465: Input Parameters:
1466: + v - the vector
1467: . n - the local size (or `PETSC_DECIDE` to have it set)
1468: - N - the global size (or `PETSC_DETERMINE` to have it set)
1470: Level: intermediate
1472: Notes:
1473: `N` cannot be `PETSC_DETERMINE` if `n` is `PETSC_DECIDE`
1475: If one processor calls this with `N` of `PETSC_DETERMINE` then all processors must, otherwise the program will hang.
1477: If `n` is not `PETSC_DECIDE`, then the value determines the `PetscLayout` of the vector and the ranges returned by
1478: `VecGetOwnershipRange()` and `VecGetOwnershipRanges()`
1480: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecCreateSeq()`, `VecCreateMPI()`, `VecGetSize()`, `PetscSplitOwnership()`, `PetscLayout`,
1481: `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`, `MatSetSizes()`
1482: @*/
1483: PetscErrorCode VecSetSizes(Vec v, PetscInt n, PetscInt N)
1484: {
1485: PetscFunctionBegin;
1487: if (N >= 0) {
1489: PetscCheck(n <= N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local size %" PetscInt_FMT " cannot be larger than global size %" PetscInt_FMT, n, N);
1490: }
1491: 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,
1492: v->map->n, v->map->N);
1493: v->map->n = n;
1494: v->map->N = N;
1495: PetscTryTypeMethod(v, create);
1496: v->ops->create = NULL;
1497: PetscFunctionReturn(PETSC_SUCCESS);
1498: }
1500: /*@
1501: VecSetBlockSize - Sets the block size for future calls to `VecSetValuesBlocked()`
1502: and `VecSetValuesBlockedLocal()`.
1504: Logically Collective
1506: Input Parameters:
1507: + v - the vector
1508: - bs - the blocksize
1510: Level: advanced
1512: Note:
1513: All vectors obtained by `VecDuplicate()` inherit the same blocksize.
1515: Vectors obtained with `DMCreateGlobalVector()` and `DMCreateLocalVector()` generally already have a blocksize set based on the state of the `DM`
1517: .seealso: [](ch_vectors), `Vec`, `VecSetValuesBlocked()`, `VecSetLocalToGlobalMapping()`, `VecGetBlockSize()`
1518: @*/
1519: PetscErrorCode VecSetBlockSize(Vec v, PetscInt bs)
1520: {
1521: PetscFunctionBegin;
1524: PetscCall(PetscLayoutSetBlockSize(v->map, bs));
1525: v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1526: PetscFunctionReturn(PETSC_SUCCESS);
1527: }
1529: /*@
1530: VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for `VecSetValuesBlocked()`
1531: and `VecSetValuesBlockedLocal()`.
1533: Not Collective
1535: Input Parameter:
1536: . v - the vector
1538: Output Parameter:
1539: . bs - the blocksize
1541: Level: advanced
1543: Note:
1544: All vectors obtained by `VecDuplicate()` inherit the same blocksize.
1546: .seealso: [](ch_vectors), `Vec`, `VecSetValuesBlocked()`, `VecSetLocalToGlobalMapping()`, `VecSetBlockSize()`
1547: @*/
1548: PetscErrorCode VecGetBlockSize(Vec v, PetscInt *bs)
1549: {
1550: PetscFunctionBegin;
1552: PetscAssertPointer(bs, 2);
1553: PetscCall(PetscLayoutGetBlockSize(v->map, bs));
1554: PetscFunctionReturn(PETSC_SUCCESS);
1555: }
1557: /*@C
1558: VecSetOptionsPrefix - Sets the prefix used for searching for all
1559: `Vec` options in the database.
1561: Logically Collective
1563: Input Parameters:
1564: + v - the `Vec` context
1565: - prefix - the prefix to prepend to all option names
1567: Level: advanced
1569: Note:
1570: A hyphen (-) must NOT be given at the beginning of the prefix name.
1571: The first character of all runtime options is AUTOMATICALLY the hyphen.
1573: .seealso: [](ch_vectors), `Vec`, `VecSetFromOptions()`
1574: @*/
1575: PetscErrorCode VecSetOptionsPrefix(Vec v, const char prefix[])
1576: {
1577: PetscFunctionBegin;
1579: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)v, prefix));
1580: PetscFunctionReturn(PETSC_SUCCESS);
1581: }
1583: /*@C
1584: VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1585: `Vec` options in the database.
1587: Logically Collective
1589: Input Parameters:
1590: + v - the `Vec` context
1591: - prefix - the prefix to prepend to all option names
1593: Level: advanced
1595: Note:
1596: A hyphen (-) must NOT be given at the beginning of the prefix name.
1597: The first character of all runtime options is AUTOMATICALLY the hyphen.
1599: .seealso: [](ch_vectors), `Vec`, `VecGetOptionsPrefix()`
1600: @*/
1601: PetscErrorCode VecAppendOptionsPrefix(Vec v, const char prefix[])
1602: {
1603: PetscFunctionBegin;
1605: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)v, prefix));
1606: PetscFunctionReturn(PETSC_SUCCESS);
1607: }
1609: /*@C
1610: VecGetOptionsPrefix - Sets the prefix used for searching for all
1611: Vec options in the database.
1613: Not Collective
1615: Input Parameter:
1616: . v - the `Vec` context
1618: Output Parameter:
1619: . prefix - pointer to the prefix string used
1621: Level: advanced
1623: Fortran Notes:
1624: The user must pass in a string `prefix` of
1625: sufficient length to hold the prefix.
1627: .seealso: [](ch_vectors), `Vec`, `VecAppendOptionsPrefix()`
1628: @*/
1629: PetscErrorCode VecGetOptionsPrefix(Vec v, const char *prefix[])
1630: {
1631: PetscFunctionBegin;
1633: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)v, prefix));
1634: PetscFunctionReturn(PETSC_SUCCESS);
1635: }
1637: /*@
1638: VecSetUp - Sets up the internal vector data structures for the later use.
1640: Collective
1642: Input Parameter:
1643: . v - the `Vec` context
1645: Level: advanced
1647: Notes:
1648: For basic use of the `Vec` classes the user need not explicitly call
1649: `VecSetUp()`, since these actions will happen automatically.
1651: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecDestroy()`
1652: @*/
1653: PetscErrorCode VecSetUp(Vec v)
1654: {
1655: PetscMPIInt size;
1657: PetscFunctionBegin;
1659: PetscCheck(v->map->n >= 0 || v->map->N >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Sizes not set");
1660: if (!((PetscObject)v)->type_name) {
1661: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1662: if (size == 1) {
1663: PetscCall(VecSetType(v, VECSEQ));
1664: } else {
1665: PetscCall(VecSetType(v, VECMPI));
1666: }
1667: }
1668: PetscFunctionReturn(PETSC_SUCCESS);
1669: }
1671: /*
1672: These currently expose the PetscScalar/PetscReal in updating the
1673: cached norm. If we push those down into the implementation these
1674: will become independent of PetscScalar/PetscReal
1675: */
1677: PetscErrorCode VecCopyAsync_Private(Vec x, Vec y, PetscDeviceContext dctx)
1678: {
1679: PetscBool flgs[4];
1680: PetscReal norms[4] = {0.0, 0.0, 0.0, 0.0};
1682: PetscFunctionBegin;
1687: if (x == y) PetscFunctionReturn(PETSC_SUCCESS);
1688: VecCheckSameLocalSize(x, 1, y, 2);
1689: VecCheckAssembled(x);
1690: PetscCall(VecSetErrorIfLocked(y, 2));
1692: #if !defined(PETSC_USE_MIXED_PRECISION)
1693: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
1694: #endif
1696: PetscCall(PetscLogEventBegin(VEC_Copy, x, y, 0, 0));
1697: #if defined(PETSC_USE_MIXED_PRECISION)
1698: extern PetscErrorCode VecGetArray(Vec, double **);
1699: extern PetscErrorCode VecRestoreArray(Vec, double **);
1700: extern PetscErrorCode VecGetArray(Vec, float **);
1701: extern PetscErrorCode VecRestoreArray(Vec, float **);
1702: extern PetscErrorCode VecGetArrayRead(Vec, const double **);
1703: extern PetscErrorCode VecRestoreArrayRead(Vec, const double **);
1704: extern PetscErrorCode VecGetArrayRead(Vec, const float **);
1705: extern PetscErrorCode VecRestoreArrayRead(Vec, const float **);
1706: if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1707: PetscInt i, n;
1708: const float *xx;
1709: double *yy;
1710: PetscCall(VecGetArrayRead(x, &xx));
1711: PetscCall(VecGetArray(y, &yy));
1712: PetscCall(VecGetLocalSize(x, &n));
1713: for (i = 0; i < n; i++) yy[i] = xx[i];
1714: PetscCall(VecRestoreArrayRead(x, &xx));
1715: PetscCall(VecRestoreArray(y, &yy));
1716: } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1717: PetscInt i, n;
1718: float *yy;
1719: const double *xx;
1720: PetscCall(VecGetArrayRead(x, &xx));
1721: PetscCall(VecGetArray(y, &yy));
1722: PetscCall(VecGetLocalSize(x, &n));
1723: for (i = 0; i < n; i++) yy[i] = (float)xx[i];
1724: PetscCall(VecRestoreArrayRead(x, &xx));
1725: PetscCall(VecRestoreArray(y, &yy));
1726: } else PetscUseTypeMethod(x, copy, y);
1727: #else
1728: VecMethodDispatch(x, dctx, VecAsyncFnName(Copy), copy, (Vec, Vec, PetscDeviceContext), y);
1729: #endif
1731: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1732: #if !defined(PETSC_USE_MIXED_PRECISION)
1733: for (PetscInt i = 0; i < 4; i++) {
1734: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)y, NormIds[i], norms[i]));
1735: }
1736: #endif
1738: PetscCall(PetscLogEventEnd(VEC_Copy, x, y, 0, 0));
1739: PetscFunctionReturn(PETSC_SUCCESS);
1740: }
1742: /*@
1743: VecCopy - Copies a vector `y = x`
1745: Logically Collective
1747: Input Parameter:
1748: . x - the vector
1750: Output Parameter:
1751: . y - the copy
1753: Level: beginner
1755: Note:
1756: For default parallel PETSc vectors, both `x` and `y` must be distributed in
1757: the same manner; local copies are done.
1759: Developer Notes:
1760: `PetscCheckSameTypeAndComm`(x,1,y,2) is not used on these vectors because we allow one
1761: of the vectors to be sequential and one to be parallel so long as both have the same
1762: local sizes. This is used in some internal functions in PETSc.
1764: .seealso: [](ch_vectors), `Vec`, `VecDuplicate()`
1765: @*/
1766: PetscErrorCode VecCopy(Vec x, Vec y)
1767: {
1768: PetscFunctionBegin;
1769: PetscCall(VecCopyAsync_Private(x, y, NULL));
1770: PetscFunctionReturn(PETSC_SUCCESS);
1771: }
1773: PetscErrorCode VecSwapAsync_Private(Vec x, Vec y, PetscDeviceContext dctx)
1774: {
1775: PetscReal normxs[4], normys[4];
1776: PetscBool flgxs[4], flgys[4];
1778: PetscFunctionBegin;
1783: PetscCheckSameTypeAndComm(x, 1, y, 2);
1784: VecCheckSameSize(x, 1, y, 2);
1785: VecCheckAssembled(x);
1786: VecCheckAssembled(y);
1787: PetscCall(VecSetErrorIfLocked(x, 1));
1788: PetscCall(VecSetErrorIfLocked(y, 2));
1790: for (PetscInt i = 0; i < 4; i++) {
1791: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], normxs[i], flgxs[i]));
1792: PetscCall(PetscObjectComposedDataGetReal((PetscObject)y, NormIds[i], normys[i], flgys[i]));
1793: }
1795: PetscCall(PetscLogEventBegin(VEC_Swap, x, y, 0, 0));
1796: VecMethodDispatch(x, dctx, VecAsyncFnName(Swap), swap, (Vec, Vec, PetscDeviceContext), y);
1797: PetscCall(PetscLogEventEnd(VEC_Swap, x, y, 0, 0));
1799: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1800: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1801: for (PetscInt i = 0; i < 4; i++) {
1802: if (flgxs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)y, NormIds[i], normxs[i]));
1803: if (flgys[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], normys[i]));
1804: }
1805: PetscFunctionReturn(PETSC_SUCCESS);
1806: }
1807: /*@
1808: VecSwap - Swaps the values between two vectors, `x` and `y`.
1810: Logically Collective
1812: Input Parameters:
1813: + x - the first vector
1814: - y - the second vector
1816: Level: advanced
1818: .seealso: [](ch_vectors), `Vec`, `VecSet()`
1819: @*/
1820: PetscErrorCode VecSwap(Vec x, Vec y)
1821: {
1822: PetscFunctionBegin;
1823: PetscCall(VecSwapAsync_Private(x, y, NULL));
1824: PetscFunctionReturn(PETSC_SUCCESS);
1825: }
1827: /*@C
1828: VecStashViewFromOptions - Processes command line options to determine if/how a `VecStash` object is to be viewed.
1830: Collective
1832: Input Parameters:
1833: + obj - the `Vec` containing a stash
1834: . bobj - optional other object that provides the prefix
1835: - optionname - option to activate viewing
1837: Level: intermediate
1839: Developer Notes:
1840: This cannot use `PetscObjectViewFromOptions()` because it takes a `Vec` as an argument but does not use `VecView()`
1842: .seealso: [](ch_vectors), `Vec`, `VecStashSetInitialSize()`
1843: @*/
1844: PetscErrorCode VecStashViewFromOptions(Vec obj, PetscObject bobj, const char optionname[])
1845: {
1846: PetscViewer viewer;
1847: PetscBool flg;
1848: PetscViewerFormat format;
1849: char *prefix;
1851: PetscFunctionBegin;
1852: prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
1853: PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)obj), ((PetscObject)obj)->options, prefix, optionname, &viewer, &format, &flg));
1854: if (flg) {
1855: PetscCall(PetscViewerPushFormat(viewer, format));
1856: PetscCall(VecStashView(obj, viewer));
1857: PetscCall(PetscViewerPopFormat(viewer));
1858: PetscCall(PetscOptionsRestoreViewer(&viewer));
1859: }
1860: PetscFunctionReturn(PETSC_SUCCESS);
1861: }
1863: /*@
1864: VecStashView - Prints the entries in the vector stash and block stash.
1866: Collective
1868: Input Parameters:
1869: + v - the vector
1870: - viewer - the viewer
1872: Level: advanced
1874: .seealso: [](ch_vectors), `Vec`, `VecSetBlockSize()`, `VecSetValues()`, `VecSetValuesBlocked()`
1875: @*/
1876: PetscErrorCode VecStashView(Vec v, PetscViewer viewer)
1877: {
1878: PetscMPIInt rank;
1879: PetscInt i, j;
1880: PetscBool match;
1881: VecStash *s;
1882: PetscScalar val;
1884: PetscFunctionBegin;
1887: PetscCheckSameComm(v, 1, viewer, 2);
1889: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &match));
1890: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_SUP, "Stash viewer only works with ASCII viewer not %s", ((PetscObject)v)->type_name);
1891: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1892: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank));
1893: s = &v->bstash;
1895: /* print block stash */
1896: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
1897: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Vector Block stash size %" PetscInt_FMT " block size %" PetscInt_FMT "\n", rank, s->n, s->bs));
1898: for (i = 0; i < s->n; i++) {
1899: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " ", rank, s->idx[i]));
1900: for (j = 0; j < s->bs; j++) {
1901: val = s->array[i * s->bs + j];
1902: #if defined(PETSC_USE_COMPLEX)
1903: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "(%18.16e %18.16e) ", (double)PetscRealPart(val), (double)PetscImaginaryPart(val)));
1904: #else
1905: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%18.16e ", (double)val));
1906: #endif
1907: }
1908: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
1909: }
1910: PetscCall(PetscViewerFlush(viewer));
1912: s = &v->stash;
1914: /* print basic stash */
1915: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Vector stash size %" PetscInt_FMT "\n", rank, s->n));
1916: for (i = 0; i < s->n; i++) {
1917: val = s->array[i];
1918: #if defined(PETSC_USE_COMPLEX)
1919: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " (%18.16e %18.16e) ", rank, s->idx[i], (double)PetscRealPart(val), (double)PetscImaginaryPart(val)));
1920: #else
1921: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " %18.16e\n", rank, s->idx[i], (double)val));
1922: #endif
1923: }
1924: PetscCall(PetscViewerFlush(viewer));
1925: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1926: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1927: PetscFunctionReturn(PETSC_SUCCESS);
1928: }
1930: PetscErrorCode PetscOptionsGetVec(PetscOptions options, const char prefix[], const char key[], Vec v, PetscBool *set)
1931: {
1932: PetscInt i, N, rstart, rend;
1933: PetscScalar *xx;
1934: PetscReal *xreal;
1935: PetscBool iset;
1937: PetscFunctionBegin;
1938: PetscCall(VecGetOwnershipRange(v, &rstart, &rend));
1939: PetscCall(VecGetSize(v, &N));
1940: PetscCall(PetscCalloc1(N, &xreal));
1941: PetscCall(PetscOptionsGetRealArray(options, prefix, key, xreal, &N, &iset));
1942: if (iset) {
1943: PetscCall(VecGetArray(v, &xx));
1944: for (i = rstart; i < rend; i++) xx[i - rstart] = xreal[i];
1945: PetscCall(VecRestoreArray(v, &xx));
1946: }
1947: PetscCall(PetscFree(xreal));
1948: if (set) *set = iset;
1949: PetscFunctionReturn(PETSC_SUCCESS);
1950: }
1952: /*@
1953: VecGetLayout - get `PetscLayout` describing a vector layout
1955: Not Collective
1957: Input Parameter:
1958: . x - the vector
1960: Output Parameter:
1961: . map - the layout
1963: Level: developer
1965: Note:
1966: The layout determines what vector elements are contained on each MPI process
1968: .seealso: [](ch_vectors), `PetscLayout`, `Vec`, `VecGetSizes()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
1969: @*/
1970: PetscErrorCode VecGetLayout(Vec x, PetscLayout *map)
1971: {
1972: PetscFunctionBegin;
1974: PetscAssertPointer(map, 2);
1975: *map = x->map;
1976: PetscFunctionReturn(PETSC_SUCCESS);
1977: }
1979: /*@
1980: VecSetLayout - set `PetscLayout` describing vector layout
1982: Not Collective
1984: Input Parameters:
1985: + x - the vector
1986: - map - the layout
1988: Level: developer
1990: Note:
1991: It is normally only valid to replace the layout with a layout known to be equivalent.
1993: .seealso: [](ch_vectors), `Vec`, `PetscLayout`, `VecGetLayout()`, `VecGetSizes()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
1994: @*/
1995: PetscErrorCode VecSetLayout(Vec x, PetscLayout map)
1996: {
1997: PetscFunctionBegin;
1999: PetscCall(PetscLayoutReference(map, &x->map));
2000: PetscFunctionReturn(PETSC_SUCCESS);
2001: }
2003: PetscErrorCode VecSetInf(Vec xin)
2004: {
2005: // use of variables one and zero over just doing 1.0/0.0 is deliberate. MSVC complains that
2006: // we are dividing by zero in the latter case (ostensibly because dividing by 0 is UB, but
2007: // only for *integers* not floats).
2008: const PetscScalar one = 1.0, zero = 0.0;
2009: PetscScalar inf;
2011: PetscFunctionBegin;
2012: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2013: inf = one / zero;
2014: PetscCall(PetscFPTrapPop());
2015: if (xin->ops->set) {
2016: PetscUseTypeMethod(xin, set, inf);
2017: } else {
2018: PetscInt n;
2019: PetscScalar *xx;
2021: PetscCall(VecGetLocalSize(xin, &n));
2022: PetscCall(VecGetArrayWrite(xin, &xx));
2023: for (PetscInt i = 0; i < n; ++i) xx[i] = inf;
2024: PetscCall(VecRestoreArrayWrite(xin, &xx));
2025: }
2026: PetscFunctionReturn(PETSC_SUCCESS);
2027: }
2029: /*@
2030: VecBindToCPU - marks a vector to temporarily stay on the CPU and perform computations on the CPU
2032: Logically collective
2034: Input Parameters:
2035: + v - the vector
2036: - flg - bind to the CPU if value of `PETSC_TRUE`
2038: Level: intermediate
2040: .seealso: [](ch_vectors), `Vec`, `VecBoundToCPU()`
2041: @*/
2042: PetscErrorCode VecBindToCPU(Vec v, PetscBool flg)
2043: {
2044: PetscFunctionBegin;
2047: #if defined(PETSC_HAVE_DEVICE)
2048: if (v->boundtocpu == flg) PetscFunctionReturn(PETSC_SUCCESS);
2049: v->boundtocpu = flg;
2050: PetscTryTypeMethod(v, bindtocpu, flg);
2051: #endif
2052: PetscFunctionReturn(PETSC_SUCCESS);
2053: }
2055: /*@
2056: VecBoundToCPU - query if a vector is bound to the CPU
2058: Not collective
2060: Input Parameter:
2061: . v - the vector
2063: Output Parameter:
2064: . flg - the logical flag
2066: Level: intermediate
2068: .seealso: [](ch_vectors), `Vec`, `VecBindToCPU()`
2069: @*/
2070: PetscErrorCode VecBoundToCPU(Vec v, PetscBool *flg)
2071: {
2072: PetscFunctionBegin;
2074: PetscAssertPointer(flg, 2);
2075: #if defined(PETSC_HAVE_DEVICE)
2076: *flg = v->boundtocpu;
2077: #else
2078: *flg = PETSC_TRUE;
2079: #endif
2080: PetscFunctionReturn(PETSC_SUCCESS);
2081: }
2083: /*@
2084: VecSetBindingPropagates - Sets whether the state of being bound to the CPU for a GPU vector type propagates to child and some other associated objects
2086: Input Parameters:
2087: + v - the vector
2088: - flg - flag indicating whether the boundtocpu flag should be propagated
2090: Level: developer
2092: Notes:
2093: 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.
2094: The created vectors will also have their bindingpropagates flag set to true.
2096: Developer Notes:
2097: If a `DMDA` has the `-dm_bind_below option` set to true, then vectors created by `DMCreateGlobalVector()` will have `VecSetBindingPropagates()` called on them to
2098: set their bindingpropagates flag to true.
2100: .seealso: [](ch_vectors), `Vec`, `MatSetBindingPropagates()`, `VecGetBindingPropagates()`
2101: @*/
2102: PetscErrorCode VecSetBindingPropagates(Vec v, PetscBool flg)
2103: {
2104: PetscFunctionBegin;
2106: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2107: v->bindingpropagates = flg;
2108: #endif
2109: PetscFunctionReturn(PETSC_SUCCESS);
2110: }
2112: /*@
2113: VecGetBindingPropagates - Gets whether the state of being bound to the CPU for a GPU vector type propagates to child and some other associated objects
2115: Input Parameter:
2116: . v - the vector
2118: Output Parameter:
2119: . flg - flag indicating whether the boundtocpu flag will be propagated
2121: Level: developer
2123: .seealso: [](ch_vectors), `Vec`, `VecSetBindingPropagates()`
2124: @*/
2125: PetscErrorCode VecGetBindingPropagates(Vec v, PetscBool *flg)
2126: {
2127: PetscFunctionBegin;
2129: PetscAssertPointer(flg, 2);
2130: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2131: *flg = v->bindingpropagates;
2132: #else
2133: *flg = PETSC_FALSE;
2134: #endif
2135: PetscFunctionReturn(PETSC_SUCCESS);
2136: }
2138: /*@C
2139: VecSetPinnedMemoryMin - Set the minimum data size for which pinned memory will be used for host (CPU) allocations.
2141: Logically Collective
2143: Input Parameters:
2144: + v - the vector
2145: - mbytes - minimum data size in bytes
2147: Options Database Key:
2148: . -vec_pinned_memory_min <size> - minimum size (in bytes) for an allocation to use pinned memory on host.
2150: Level: developer
2152: Note:
2153: Specifying -1 ensures that pinned memory will never be used.
2155: .seealso: [](ch_vectors), `Vec`, `VecGetPinnedMemoryMin()`
2156: @*/
2157: PetscErrorCode VecSetPinnedMemoryMin(Vec v, size_t mbytes)
2158: {
2159: PetscFunctionBegin;
2161: #if PetscDefined(HAVE_DEVICE)
2162: v->minimum_bytes_pinned_memory = mbytes;
2163: #endif
2164: PetscFunctionReturn(PETSC_SUCCESS);
2165: }
2167: /*@C
2168: VecGetPinnedMemoryMin - Get the minimum data size for which pinned memory will be used for host (CPU) allocations.
2170: Logically Collective
2172: Input Parameter:
2173: . v - the vector
2175: Output Parameter:
2176: . mbytes - minimum data size in bytes
2178: Level: developer
2180: .seealso: [](ch_vectors), `Vec`, `VecSetPinnedMemoryMin()`
2181: @*/
2182: PetscErrorCode VecGetPinnedMemoryMin(Vec v, size_t *mbytes)
2183: {
2184: PetscFunctionBegin;
2186: PetscAssertPointer(mbytes, 2);
2187: #if PetscDefined(HAVE_DEVICE)
2188: *mbytes = v->minimum_bytes_pinned_memory;
2189: #endif
2190: PetscFunctionReturn(PETSC_SUCCESS);
2191: }
2193: /*@
2194: VecGetOffloadMask - Get the offload mask of a `Vec`
2196: Not Collective
2198: Input Parameter:
2199: . v - the vector
2201: Output Parameter:
2202: . mask - corresponding `PetscOffloadMask` enum value.
2204: Level: intermediate
2206: .seealso: [](ch_vectors), `Vec`, `VecCreateSeqCUDA()`, `VecCreateSeqViennaCL()`, `VecGetArray()`, `VecGetType()`
2207: @*/
2208: PetscErrorCode VecGetOffloadMask(Vec v, PetscOffloadMask *mask)
2209: {
2210: PetscFunctionBegin;
2212: PetscAssertPointer(mask, 2);
2213: *mask = v->offloadmask;
2214: PetscFunctionReturn(PETSC_SUCCESS);
2215: }
2217: #if !defined(PETSC_HAVE_VIENNACL)
2218: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLContext(Vec v, PETSC_UINTPTR_T *ctx)
2219: {
2220: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_context");
2221: }
2223: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLQueue(Vec v, PETSC_UINTPTR_T *queue)
2224: {
2225: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_command_queue");
2226: }
2228: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMem(Vec v, PETSC_UINTPTR_T *queue)
2229: {
2230: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2231: }
2233: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMemRead(Vec v, PETSC_UINTPTR_T *queue)
2234: {
2235: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2236: }
2238: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMemWrite(Vec v, PETSC_UINTPTR_T *queue)
2239: {
2240: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2241: }
2243: PETSC_EXTERN PetscErrorCode VecViennaCLRestoreCLMemWrite(Vec v)
2244: {
2245: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to restore a Vec's cl_mem");
2246: }
2247: #endif
2249: 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)
2250: {
2251: const PetscScalar *u, *y;
2252: const PetscScalar *atola = NULL, *rtola = NULL, *erra = NULL;
2253: PetscInt n, n_loc = 0, na_loc = 0, nr_loc = 0;
2254: PetscReal nrm = 0, nrma = 0, nrmr = 0, err_loc[6];
2256: PetscFunctionBegin;
2257: #define SkipSmallValue(a, b, tol) \
2258: if (PetscAbsScalar(a) < tol || PetscAbsScalar(b) < tol) continue
2260: PetscCall(VecGetLocalSize(U, &n));
2261: PetscCall(VecGetArrayRead(U, &u));
2262: PetscCall(VecGetArrayRead(Y, &y));
2263: if (E) PetscCall(VecGetArrayRead(E, &erra));
2264: if (vatol) PetscCall(VecGetArrayRead(vatol, &atola));
2265: if (vrtol) PetscCall(VecGetArrayRead(vrtol, &rtola));
2266: for (PetscInt i = 0; i < n; i++) {
2267: PetscReal err, tol, tola, tolr;
2269: SkipSmallValue(y[i], u[i], ignore_max);
2270: atol = atola ? PetscRealPart(atola[i]) : atol;
2271: rtol = rtola ? PetscRealPart(rtola[i]) : rtol;
2272: err = erra ? PetscAbsScalar(erra[i]) : PetscAbsScalar(y[i] - u[i]);
2273: tola = atol;
2274: tolr = rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
2275: tol = tola + tolr;
2276: if (tola > 0.) {
2277: if (wnormtype == NORM_INFINITY) nrma = PetscMax(nrma, err / tola);
2278: else nrma += PetscSqr(err / tola);
2279: na_loc++;
2280: }
2281: if (tolr > 0.) {
2282: if (wnormtype == NORM_INFINITY) nrmr = PetscMax(nrmr, err / tolr);
2283: else nrmr += PetscSqr(err / tolr);
2284: nr_loc++;
2285: }
2286: if (tol > 0.) {
2287: if (wnormtype == NORM_INFINITY) nrm = PetscMax(nrm, err / tol);
2288: else nrm += PetscSqr(err / tol);
2289: n_loc++;
2290: }
2291: }
2292: if (E) PetscCall(VecRestoreArrayRead(E, &erra));
2293: if (vatol) PetscCall(VecRestoreArrayRead(vatol, &atola));
2294: if (vrtol) PetscCall(VecRestoreArrayRead(vrtol, &rtola));
2295: PetscCall(VecRestoreArrayRead(U, &u));
2296: PetscCall(VecRestoreArrayRead(Y, &y));
2297: #undef SkipSmallValue
2299: err_loc[0] = nrm;
2300: err_loc[1] = nrma;
2301: err_loc[2] = nrmr;
2302: err_loc[3] = (PetscReal)n_loc;
2303: err_loc[4] = (PetscReal)na_loc;
2304: err_loc[5] = (PetscReal)nr_loc;
2305: if (wnormtype == NORM_2) {
2306: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, err_loc, 6, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)U)));
2307: } else {
2308: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, err_loc, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)U)));
2309: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, err_loc + 3, 3, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)U)));
2310: }
2311: if (wnormtype == NORM_2) {
2312: *norm = PetscSqrtReal(err_loc[0]);
2313: *norma = PetscSqrtReal(err_loc[1]);
2314: *normr = PetscSqrtReal(err_loc[2]);
2315: } else {
2316: *norm = err_loc[0];
2317: *norma = err_loc[1];
2318: *normr = err_loc[2];
2319: }
2320: *norm_loc = (PetscInt)err_loc[3];
2321: *norma_loc = (PetscInt)err_loc[4];
2322: *normr_loc = (PetscInt)err_loc[5];
2323: PetscFunctionReturn(PETSC_SUCCESS);
2324: }
2326: /*@
2327: VecErrorWeightedNorms - compute a weighted norm of the difference between two vectors
2329: Collective
2331: Input Parameters:
2332: + U - first vector to be compared
2333: . Y - second vector to be compared
2334: . E - optional third vector representing the error (if not provided, the error is ||U-Y||)
2335: . wnormtype - norm type
2336: . atol - scalar for absolute tolerance
2337: . vatol - vector representing per-entry absolute tolerances (can be ``NULL``)
2338: . rtol - scalar for relative tolerance
2339: . vrtol - vector representing per-entry relative tolerances (can be ``NULL``)
2340: - ignore_max - ignore values smaller then this value in absolute terms.
2342: Output Parameters:
2343: + norm - weighted norm
2344: . norm_loc - number of vector locations used for the weighted norm
2345: . norma - weighted norm based on the absolute tolerance
2346: . norma_loc - number of vector locations used for the absolute weighted norm
2347: . normr - weighted norm based on the relative tolerance
2348: - normr_loc - number of vector locations used for the relative weighted norm
2350: Level: developer
2352: Notes:
2353: This is primarily used for computing weighted local truncation errors in ``TS``.
2355: .seealso: [](ch_vectors), `Vec`, `NormType`, ``TSErrorWeightedNorm()``, ``TSErrorWeightedENorm()``
2356: @*/
2357: 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)
2358: {
2359: PetscFunctionBegin;
2364: if (E) {
2367: }
2370: if (vatol) {
2373: }
2375: if (vrtol) {
2378: }
2380: PetscAssertPointer(norm, 10);
2381: PetscAssertPointer(norm_loc, 11);
2382: PetscAssertPointer(norma, 12);
2383: PetscAssertPointer(norma_loc, 13);
2384: PetscAssertPointer(normr, 14);
2385: PetscAssertPointer(normr_loc, 15);
2386: PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)U), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
2388: /* There are potentially 5 vectors involved, some of them may happen to be of different type or bound to cpu.
2389: Here we check that they all implement the same operation and call it if so.
2390: Otherwise, we call the _Basic implementation that always works (provided VecGetArrayRead is implemented). */
2391: PetscBool sameop = (PetscBool)(U->ops->errorwnorm && U->ops->errorwnorm == Y->ops->errorwnorm);
2392: if (sameop && E) sameop = (PetscBool)(U->ops->errorwnorm == E->ops->errorwnorm);
2393: if (sameop && vatol) sameop = (PetscBool)(U->ops->errorwnorm == vatol->ops->errorwnorm);
2394: if (sameop && vrtol) sameop = (PetscBool)(U->ops->errorwnorm == vrtol->ops->errorwnorm);
2395: if (sameop) PetscUseTypeMethod(U, errorwnorm, Y, E, wnormtype, atol, vatol, rtol, vrtol, ignore_max, norm, norm_loc, norma, norma_loc, normr, normr_loc);
2396: else PetscCall(VecErrorWeightedNorms_Basic(U, Y, E, wnormtype, atol, vatol, rtol, vrtol, ignore_max, norm, norm_loc, norma, norma_loc, normr, normr_loc));
2397: PetscFunctionReturn(PETSC_SUCCESS);
2398: }