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) {
335: PetscCall((*async_fn)(w, x, y, dctx));
336: } else {
337: PetscCall((*pointwise_op)(w, x, y));
338: }
339: if (event) PetscCall(PetscLogEventEnd(event, x, y, w, 0));
340: PetscCall(PetscObjectStateIncrease((PetscObject)w));
341: PetscFunctionReturn(PETSC_SUCCESS);
342: }
344: PetscErrorCode VecPointwiseMaxAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
345: {
346: PetscFunctionBegin;
347: // REVIEW ME: no log event?
348: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMax), w->ops->pointwisemax));
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: /*@
353: VecPointwiseMax - Computes the component-wise maximum `w[i] = max(x[i], y[i])`.
355: Logically Collective
357: Input Parameters:
358: + x - the first input vector
359: - y - the second input vector
361: Output Parameter:
362: . w - the result
364: Level: advanced
366: Notes:
367: Any subset of the `x`, `y`, and `w` may be the same vector.
369: For complex numbers compares only the real part
371: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
372: @*/
373: PetscErrorCode VecPointwiseMax(Vec w, Vec x, Vec y)
374: {
375: PetscFunctionBegin;
376: PetscCall(VecPointwiseMaxAsync_Private(w, x, y, NULL));
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: PetscErrorCode VecPointwiseMinAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
381: {
382: PetscFunctionBegin;
383: // REVIEW ME: no log event?
384: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMin), w->ops->pointwisemin));
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: /*@
389: VecPointwiseMin - Computes the component-wise minimum `w[i] = min(x[i], y[i])`.
391: Logically Collective
393: Input Parameters:
394: + x - the first input vector
395: - y - the second input vector
397: Output Parameter:
398: . w - the result
400: Level: advanced
402: Notes:
403: Any subset of the `x`, `y`, and `w` may be the same vector.
405: For complex numbers compares only the real part
407: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
408: @*/
409: PetscErrorCode VecPointwiseMin(Vec w, Vec x, Vec y)
410: {
411: PetscFunctionBegin;
412: PetscCall(VecPointwiseMinAsync_Private(w, x, y, NULL));
413: PetscFunctionReturn(PETSC_SUCCESS);
414: }
416: PetscErrorCode VecPointwiseMaxAbsAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
417: {
418: PetscFunctionBegin;
419: // REVIEW ME: no log event?
420: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMaxAbs), w->ops->pointwisemaxabs));
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: /*@
425: VecPointwiseMaxAbs - Computes the component-wise maximum of the absolute values `w[i] = max(abs(x[i]), abs(y[i]))`.
427: Logically Collective
429: Input Parameters:
430: + x - the first input vector
431: - y - the second input vector
433: Output Parameter:
434: . w - the result
436: Level: advanced
438: Notes:
439: Any subset of the `x`, `y`, and `w` may be the same vector.
441: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMin()`, `VecPointwiseMax()`, `VecMaxPointwiseDivide()`
442: @*/
443: PetscErrorCode VecPointwiseMaxAbs(Vec w, Vec x, Vec y)
444: {
445: PetscFunctionBegin;
446: PetscCall(VecPointwiseMaxAbsAsync_Private(w, x, y, NULL));
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: PetscErrorCode VecPointwiseDivideAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
451: {
452: PetscFunctionBegin;
453: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, VEC_PointwiseDivide, VecAsyncFnName(PointwiseDivide), w->ops->pointwisedivide));
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: /*@
458: VecPointwiseDivide - Computes the component-wise division `w[i] = x[i] / y[i]`.
460: Logically Collective
462: Input Parameters:
463: + x - the numerator vector
464: - y - the denominator vector
466: Output Parameter:
467: . w - the result
469: Level: advanced
471: Note:
472: Any subset of the `x`, `y`, and `w` may be the same vector.
474: .seealso: [](ch_vectors), `Vec`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
475: @*/
476: PetscErrorCode VecPointwiseDivide(Vec w, Vec x, Vec y)
477: {
478: PetscFunctionBegin;
479: PetscCall(VecPointwiseDivideAsync_Private(w, x, y, NULL));
480: PetscFunctionReturn(PETSC_SUCCESS);
481: }
483: PetscErrorCode VecPointwiseMultAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
484: {
485: PetscFunctionBegin;
487: PetscCall(VecPointwiseApply_Private(w, x, y, dctx, VEC_PointwiseMult, VecAsyncFnName(PointwiseMult), w->ops->pointwisemult));
488: PetscFunctionReturn(PETSC_SUCCESS);
489: }
491: /*@
492: VecPointwiseMult - Computes the component-wise multiplication `w[i] = x[i] * y[i]`.
494: Logically Collective
496: Input Parameters:
497: + x - the first vector
498: - y - the second vector
500: Output Parameter:
501: . w - the result
503: Level: advanced
505: Note:
506: Any subset of the `x`, `y`, and `w` may be the same vector.
508: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
509: @*/
510: PetscErrorCode VecPointwiseMult(Vec w, Vec x, Vec y)
511: {
512: PetscFunctionBegin;
513: PetscCall(VecPointwiseMultAsync_Private(w, x, y, NULL));
514: PetscFunctionReturn(PETSC_SUCCESS);
515: }
517: /*@
518: VecDuplicate - Creates a new vector of the same type as an existing vector.
520: Collective
522: Input Parameter:
523: . v - a vector to mimic
525: Output Parameter:
526: . newv - location to put new vector
528: Level: beginner
530: Notes:
531: `VecDuplicate()` DOES NOT COPY the vector entries, but rather allocates storage
532: for the new vector. Use `VecCopy()` to copy a vector.
534: Use `VecDestroy()` to free the space. Use `VecDuplicateVecs()` to get several
535: vectors.
537: .seealso: [](ch_vectors), `Vec`, `VecDestroy()`, `VecDuplicateVecs()`, `VecCreate()`, `VecCopy()`
538: @*/
539: PetscErrorCode VecDuplicate(Vec v, Vec *newv)
540: {
541: PetscFunctionBegin;
543: PetscAssertPointer(newv, 2);
545: PetscUseTypeMethod(v, duplicate, newv);
546: #if PetscDefined(HAVE_DEVICE)
547: if (v->boundtocpu && v->bindingpropagates) {
548: PetscCall(VecSetBindingPropagates(*newv, PETSC_TRUE));
549: PetscCall(VecBindToCPU(*newv, PETSC_TRUE));
550: }
551: #endif
552: PetscCall(PetscObjectStateIncrease((PetscObject)*newv));
553: PetscFunctionReturn(PETSC_SUCCESS);
554: }
556: /*@
557: VecDestroy - Destroys a vector.
559: Collective
561: Input Parameter:
562: . v - the vector
564: Level: beginner
566: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecDuplicate()`, `VecDestroyVecs()`
567: @*/
568: PetscErrorCode VecDestroy(Vec *v)
569: {
570: PetscFunctionBegin;
571: PetscAssertPointer(v, 1);
572: if (!*v) PetscFunctionReturn(PETSC_SUCCESS);
574: if (--((PetscObject)*v)->refct > 0) {
575: *v = NULL;
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }
579: PetscCall(PetscObjectSAWsViewOff((PetscObject)*v));
580: /* destroy the internal part */
581: PetscTryTypeMethod(*v, destroy);
582: PetscCall(PetscFree((*v)->defaultrandtype));
583: /* destroy the external/common part */
584: PetscCall(PetscLayoutDestroy(&(*v)->map));
585: PetscCall(PetscHeaderDestroy(v));
586: PetscFunctionReturn(PETSC_SUCCESS);
587: }
589: /*@C
590: VecDuplicateVecs - Creates several vectors of the same type as an existing vector.
592: Collective
594: Input Parameters:
595: + m - the number of vectors to obtain
596: - v - a vector to mimic
598: Output Parameter:
599: . V - location to put pointer to array of vectors
601: Level: intermediate
603: Notes:
604: Use `VecDestroyVecs()` to free the space. Use `VecDuplicate()` to form a single
605: vector.
607: Some implementations ensure that the arrays accessed by each vector are contiguous in memory. Certain `VecMDot()` and `VecMAXPY()`
608: implementations utilize this property to use BLAS 2 operations for higher efficiency. This is especially useful in `KSPGMRES`, see
609: `KSPGMRESSetPreAllocateVectors()`.
611: Fortran Note:
612: .vb
613: Vec, pointer :: V(:)
614: .ve
616: .seealso: [](ch_vectors), `Vec`, [](ch_fortran), `VecDestroyVecs()`, `VecDuplicate()`, `VecCreate()`, `VecMDot()`, `VecMAXPY()`, `KSPGMRES`,
617: `KSPGMRESSetPreAllocateVectors()`
618: @*/
619: PetscErrorCode VecDuplicateVecs(Vec v, PetscInt m, Vec *V[])
620: {
621: PetscFunctionBegin;
623: PetscAssertPointer(V, 3);
625: PetscUseTypeMethod(v, duplicatevecs, m, V);
626: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
627: if (v->boundtocpu && v->bindingpropagates) {
628: PetscInt i;
630: for (i = 0; i < m; i++) {
631: /* Since ops->duplicatevecs might itself propagate the value of boundtocpu,
632: * avoid unnecessary overhead by only calling VecBindToCPU() if the vector isn't already bound. */
633: if (!(*V)[i]->boundtocpu) {
634: PetscCall(VecSetBindingPropagates((*V)[i], PETSC_TRUE));
635: PetscCall(VecBindToCPU((*V)[i], PETSC_TRUE));
636: }
637: }
638: }
639: #endif
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: /*@C
644: VecDestroyVecs - Frees a block of vectors obtained with `VecDuplicateVecs()`.
646: Collective
648: Input Parameters:
649: + m - the number of vectors previously obtained, if zero no vectors are destroyed
650: - vv - pointer to pointer to array of vector pointers, if `NULL` no vectors are destroyed
652: Level: intermediate
654: .seealso: [](ch_vectors), `Vec`, [](ch_fortran), `VecDuplicateVecs()`, `VecDestroyVecsf90()`
655: @*/
656: PetscErrorCode VecDestroyVecs(PetscInt m, Vec *vv[])
657: {
658: PetscFunctionBegin;
659: PetscAssertPointer(vv, 2);
660: PetscCheck(m >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Trying to destroy negative number of vectors %" PetscInt_FMT, m);
661: if (!m || !*vv) {
662: *vv = NULL;
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
667: PetscCall((*(**vv)->ops->destroyvecs)(m, *vv));
668: *vv = NULL;
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: /*@
673: VecViewFromOptions - View a vector based on values in the options database
675: Collective
677: Input Parameters:
678: + A - the vector
679: . obj - optional object that provides the options prefix for this viewing, use 'NULL' to use the prefix of `A`
680: - name - command line option
682: Level: intermediate
684: Note:
685: See `PetscObjectViewFromOptions()` to see the `PetscViewer` and PetscViewerFormat` available
687: .seealso: [](ch_vectors), `Vec`, `VecView`, `PetscObjectViewFromOptions()`, `VecCreate()`
688: @*/
689: PetscErrorCode VecViewFromOptions(Vec A, PeOp PetscObject obj, const char name[])
690: {
691: PetscFunctionBegin;
693: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
694: PetscFunctionReturn(PETSC_SUCCESS);
695: }
697: /*@
698: VecView - Views a vector object.
700: Collective
702: Input Parameters:
703: + vec - the vector
704: - viewer - an optional `PetscViewer` visualization context
706: Level: beginner
708: Notes:
709: The available visualization contexts include
710: + `PETSC_VIEWER_STDOUT_SELF` - for sequential vectors
711: . `PETSC_VIEWER_STDOUT_WORLD` - for parallel vectors created on `PETSC_COMM_WORLD`
712: - `PETSC_VIEWER_STDOUT`_(comm) - for parallel vectors created on MPI communicator comm
714: You can change the format the vector is printed using the
715: option `PetscViewerPushFormat()`.
717: The user can open alternative viewers with
718: + `PetscViewerASCIIOpen()` - Outputs vector to a specified file
719: . `PetscViewerBinaryOpen()` - Outputs vector in binary to a
720: specified file; corresponding input uses `VecLoad()`
721: . `PetscViewerDrawOpen()` - Outputs vector to an X window display
722: . `PetscViewerSocketOpen()` - Outputs vector to Socket viewer
723: - `PetscViewerHDF5Open()` - Outputs vector to HDF5 file viewer
725: The user can call `PetscViewerPushFormat()` to specify the output
726: format of ASCII printed objects (when using `PETSC_VIEWER_STDOUT_SELF`,
727: `PETSC_VIEWER_STDOUT_WORLD` and `PetscViewerASCIIOpen()`). Available formats include
728: + `PETSC_VIEWER_DEFAULT` - default, prints vector contents
729: . `PETSC_VIEWER_ASCII_MATLAB` - prints vector contents in MATLAB format
730: . `PETSC_VIEWER_ASCII_INDEX` - prints vector contents, including indices of vector elements
731: - `PETSC_VIEWER_ASCII_COMMON` - prints vector contents, using a
732: format common among all vector types
734: You can pass any number of vector objects, or other PETSc objects to the same viewer.
736: In the debugger you can do call `VecView`(v,0) to display the vector. (The same holds for any PETSc object viewer).
738: Notes for binary viewer:
739: If you pass multiple vectors to a binary viewer you can read them back in the same order
740: with `VecLoad()`.
742: If the blocksize of the vector is greater than one then you must provide a unique prefix to
743: the vector with `PetscObjectSetOptionsPrefix`((`PetscObject`)vec,"uniqueprefix"); BEFORE calling `VecView()` on the
744: vector to be stored and then set that same unique prefix on the vector that you pass to `VecLoad()`. The blocksize
745: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
746: filename. If you copy the binary file, make sure you copy the associated .info file with it.
748: See the manual page for `VecLoad()` on the exact format the binary viewer stores
749: the values in the file.
751: Notes for HDF5 Viewer:
752: The name of the `Vec` (given with `PetscObjectSetName()` is the name that is used
753: for the object in the HDF5 file. If you wish to store the same Vec into multiple
754: datasets in the same file (typically with different values), you must change its
755: name each time before calling the `VecView()`. To load the same vector,
756: the name of the Vec object passed to `VecLoad()` must be the same.
758: If the block size of the vector is greater than 1 then it is used as the first dimension in the HDF5 array.
759: If the function `PetscViewerHDF5SetBaseDimension2()`is called then even if the block size is one it will
760: be used as the first dimension in the HDF5 array (that is the HDF5 array will always be two dimensional)
761: See also `PetscViewerHDF5SetTimestep()` which adds an additional complication to reading and writing `Vec`
762: with the HDF5 viewer.
764: .seealso: [](ch_vectors), `Vec`, `VecViewFromOptions()`, `PetscViewerASCIIOpen()`, `PetscViewerDrawOpen()`, `PetscDrawLGCreate()`,
765: `PetscViewerSocketOpen()`, `PetscViewerBinaryOpen()`, `VecLoad()`, `PetscViewerCreate()`,
766: `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`, `PetscViewerHDF5SetTimestep()`
767: @*/
768: PetscErrorCode VecView(Vec vec, PetscViewer viewer)
769: {
770: PetscBool isascii;
771: PetscViewerFormat format;
772: PetscMPIInt size;
774: PetscFunctionBegin;
777: VecCheckAssembled(vec);
778: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec), &viewer));
780: PetscCall(PetscViewerGetFormat(viewer, &format));
781: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
782: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) PetscFunctionReturn(PETSC_SUCCESS);
784: PetscCheck(!vec->stash.n && !vec->bstash.n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call VecAssemblyBegin/End() before viewing this vector");
786: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
787: if (isascii) {
788: PetscInt rows, bs;
790: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)vec, viewer));
791: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
792: PetscCall(PetscViewerASCIIPushTab(viewer));
793: PetscCall(VecGetSize(vec, &rows));
794: PetscCall(VecGetBlockSize(vec, &bs));
795: if (bs != 1) {
796: PetscCall(PetscViewerASCIIPrintf(viewer, "length=%" PetscInt_FMT ", bs=%" PetscInt_FMT "\n", rows, bs));
797: } else {
798: PetscCall(PetscViewerASCIIPrintf(viewer, "length=%" PetscInt_FMT "\n", rows));
799: }
800: PetscCall(PetscViewerASCIIPopTab(viewer));
801: }
802: }
803: PetscCall(VecLockReadPush(vec));
804: PetscCall(PetscLogEventBegin(VEC_View, vec, viewer, 0, 0));
805: if ((format == PETSC_VIEWER_NATIVE || format == PETSC_VIEWER_LOAD_BALANCE) && vec->ops->viewnative) {
806: PetscUseTypeMethod(vec, viewnative, viewer);
807: } else {
808: PetscUseTypeMethod(vec, view, viewer);
809: }
810: PetscCall(VecLockReadPop(vec));
811: PetscCall(PetscLogEventEnd(VEC_View, vec, viewer, 0, 0));
812: PetscFunctionReturn(PETSC_SUCCESS);
813: }
815: #if defined(PETSC_USE_DEBUG)
816: #include <../src/sys/totalview/tv_data_display.h>
817: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
818: {
819: const PetscScalar *values;
820: char type[32];
822: TV_add_row("Local rows", "int", &v->map->n);
823: TV_add_row("Global rows", "int", &v->map->N);
824: TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
825: PetscCall(VecGetArrayRead((Vec)v, &values));
826: PetscCall(PetscSNPrintf(type, 32, "double[%" PetscInt_FMT "]", v->map->n));
827: TV_add_row("values", type, values);
828: PetscCall(VecRestoreArrayRead((Vec)v, &values));
829: return TV_format_OK;
830: }
831: #endif
833: /*@C
834: VecViewNative - Views a vector object with the original type specific viewer
836: Collective
838: Input Parameters:
839: + vec - the vector
840: - viewer - an optional `PetscViewer` visualization context
842: Level: developer
844: Note:
845: This can be used with, for example, vectors obtained with `DMCreateGlobalVector()` for a `DMDA` to display the vector
846: in the PETSc storage format (each MPI process values follow the previous MPI processes) instead of the "natural" grid
847: ordering.
849: .seealso: [](ch_vectors), `Vec`, `PetscViewerASCIIOpen()`, `PetscViewerDrawOpen()`, `PetscDrawLGCreate()`, `VecView()`
850: `PetscViewerSocketOpen()`, `PetscViewerBinaryOpen()`, `VecLoad()`, `PetscViewerCreate()`,
851: `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`, `PetscViewerHDF5SetTimestep()`
852: @*/
853: PetscErrorCode VecViewNative(Vec vec, PetscViewer viewer)
854: {
855: PetscFunctionBegin;
858: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec), &viewer));
860: PetscUseTypeMethod(vec, viewnative, viewer);
861: PetscFunctionReturn(PETSC_SUCCESS);
862: }
864: /*@
865: VecGetSize - Returns the global number of elements of the vector.
867: Not Collective
869: Input Parameter:
870: . x - the vector
872: Output Parameter:
873: . size - the global length of the vector
875: Level: beginner
877: .seealso: [](ch_vectors), `Vec`, `VecGetLocalSize()`
878: @*/
879: PetscErrorCode VecGetSize(Vec x, PetscInt *size)
880: {
881: PetscFunctionBegin;
883: PetscAssertPointer(size, 2);
885: PetscUseTypeMethod(x, getsize, size);
886: PetscFunctionReturn(PETSC_SUCCESS);
887: }
889: /*@
890: VecGetLocalSize - Returns the number of elements of the vector stored
891: in local memory (that is on this MPI process)
893: Not Collective
895: Input Parameter:
896: . x - the vector
898: Output Parameter:
899: . size - the length of the local piece of the vector
901: Level: beginner
903: .seealso: [](ch_vectors), `Vec`, `VecGetSize()`
904: @*/
905: PetscErrorCode VecGetLocalSize(Vec x, PetscInt *size)
906: {
907: PetscFunctionBegin;
909: PetscAssertPointer(size, 2);
911: PetscUseTypeMethod(x, getlocalsize, size);
912: PetscFunctionReturn(PETSC_SUCCESS);
913: }
915: /*@
916: VecGetOwnershipRange - Returns the range of indices owned by
917: this process. The vector is laid out with the
918: first `n1` elements on the first processor, next `n2` elements on the
919: second, etc. For certain parallel layouts this range may not be
920: well defined.
922: Not Collective
924: Input Parameter:
925: . x - the vector
927: Output Parameters:
928: + low - the first local element, pass in `NULL` if not interested
929: - high - one more than the last local element, pass in `NULL` if not interested
931: Level: beginner
933: Notes:
934: If the `Vec` was obtained from a `DM` with `DMCreateGlobalVector()`, then the range values are determined by the specific `DM`.
936: If the `Vec` was created directly the range values are determined by the local size passed to `VecSetSizes()` or `VecCreateMPI()`.
937: If `PETSC_DECIDE` was passed as the local size, then the vector uses default values for the range using `PetscSplitOwnership()`.
939: The high argument is one more than the last element stored locally.
941: For certain `DM`, such as `DMDA`, it is better to use `DM` specific routines, such as `DMDAGetGhostCorners()`, to determine
942: the local values in the vector.
944: .seealso: [](ch_vectors), `Vec`, `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `VecGetOwnershipRanges()`, `PetscSplitOwnership()`,
945: `VecSetSizes()`, `VecCreateMPI()`, `PetscLayout`, `DMDAGetGhostCorners()`, `DM`
946: @*/
947: PetscErrorCode VecGetOwnershipRange(Vec x, PetscInt *low, PetscInt *high)
948: {
949: PetscFunctionBegin;
952: if (low) PetscAssertPointer(low, 2);
953: if (high) PetscAssertPointer(high, 3);
954: if (low) *low = x->map->rstart;
955: if (high) *high = x->map->rend;
956: PetscFunctionReturn(PETSC_SUCCESS);
957: }
959: /*@C
960: VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
961: The vector is laid out with the
962: first `n1` elements on the first processor, next `n2` elements on the
963: second, etc. For certain parallel layouts this range may not be
964: well defined.
966: Not Collective
968: Input Parameter:
969: . x - the vector
971: Output Parameter:
972: . ranges - array of length `size` + 1 with the start and end+1 for each process
974: Level: beginner
976: Notes:
977: If the `Vec` was obtained from a `DM` with `DMCreateGlobalVector()`, then the range values are determined by the specific `DM`.
979: If the `Vec` was created directly the range values are determined by the local size passed to `VecSetSizes()` or `VecCreateMPI()`.
980: If `PETSC_DECIDE` was passed as the local size, then the vector uses default values for the range using `PetscSplitOwnership()`.
982: The high argument is one more than the last element stored locally.
984: For certain `DM`, such as `DMDA`, it is better to use `DM` specific routines, such as `DMDAGetGhostCorners()`, to determine
985: the local values in the vector.
987: The high argument is one more than the last element stored locally.
989: If `ranges` are used after all vectors that share the ranges has been destroyed, then the program will crash accessing `ranges`.
991: Fortran Note:
992: The argument `ranges` must be declared as
993: .vb
994: PetscInt, pointer :: ranges(:)
995: .ve
996: and you have to return it with a call to `VecRestoreOwnershipRanges()` when no longer needed
998: .seealso: [](ch_vectors), `Vec`, `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `VecGetOwnershipRange()`, `PetscSplitOwnership()`,
999: `VecSetSizes()`, `VecCreateMPI()`, `PetscLayout`, `DMDAGetGhostCorners()`, `DM`
1000: @*/
1001: PetscErrorCode VecGetOwnershipRanges(Vec x, const PetscInt *ranges[])
1002: {
1003: PetscFunctionBegin;
1006: PetscCall(PetscLayoutGetRanges(x->map, ranges));
1007: PetscFunctionReturn(PETSC_SUCCESS);
1008: }
1010: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
1011: /*@
1012: VecSetOption - Sets an option for controlling a vector's behavior.
1014: Collective
1016: Input Parameters:
1017: + x - the vector
1018: . op - the option
1019: - flag - turn the option on or off
1021: Supported Options:
1022: + `VEC_IGNORE_OFF_PROC_ENTRIES` - which causes `VecSetValues()` to ignore
1023: entries destined to be stored on a separate processor. This can be used
1024: to eliminate the global reduction in the `VecAssemblyBegin()` if you know
1025: that you have only used `VecSetValues()` to set local elements
1026: . `VEC_IGNORE_NEGATIVE_INDICES` - which means you can pass negative indices
1027: in ix in calls to `VecSetValues()` or `VecGetValues()`. These rows are simply
1028: ignored.
1029: - `VEC_SUBSET_OFF_PROC_ENTRIES` - which causes `VecAssemblyBegin()` to assume that the off-process
1030: entries will always be a subset (possibly equal) of the off-process entries set on the
1031: first assembly which had a true `VEC_SUBSET_OFF_PROC_ENTRIES` and the vector has not
1032: changed this flag afterwards. If this assembly is not such first assembly, then this
1033: assembly can reuse the communication pattern setup in that first assembly, thus avoiding
1034: a global reduction. Subsequent assemblies setting off-process values should use the same
1035: InsertMode as the first assembly.
1037: Level: intermediate
1039: Developer Notes:
1040: The `InsertMode` restriction could be removed by packing the stash messages out of place.
1042: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`
1043: @*/
1044: PetscErrorCode VecSetOption(Vec x, VecOption op, PetscBool flag)
1045: {
1046: PetscFunctionBegin;
1049: PetscTryTypeMethod(x, setoption, op, flag);
1050: PetscFunctionReturn(PETSC_SUCCESS);
1051: }
1053: /* Default routines for obtaining and releasing; */
1054: /* may be used by any implementation */
1055: PetscErrorCode VecDuplicateVecs_Default(Vec w, PetscInt m, Vec *V[])
1056: {
1057: PetscFunctionBegin;
1058: PetscCheck(m > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "m must be > 0: m = %" PetscInt_FMT, m);
1059: PetscCall(PetscMalloc1(m, V));
1060: for (PetscInt i = 0; i < m; i++) PetscCall(VecDuplicate(w, *V + i));
1061: PetscFunctionReturn(PETSC_SUCCESS);
1062: }
1064: PetscErrorCode VecDestroyVecs_Default(PetscInt m, Vec v[])
1065: {
1066: PetscInt i;
1068: PetscFunctionBegin;
1069: PetscAssertPointer(v, 2);
1070: for (i = 0; i < m; i++) PetscCall(VecDestroy(&v[i]));
1071: PetscCall(PetscFree(v));
1072: PetscFunctionReturn(PETSC_SUCCESS);
1073: }
1075: /*@
1076: VecResetArray - Resets a vector to use its default memory. Call this
1077: after the use of `VecPlaceArray()`.
1079: Not Collective
1081: Input Parameter:
1082: . vec - the vector
1084: Level: developer
1086: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecPlaceArray()`
1087: @*/
1088: PetscErrorCode VecResetArray(Vec vec)
1089: {
1090: PetscFunctionBegin;
1093: PetscUseTypeMethod(vec, resetarray);
1094: PetscCall(PetscObjectStateIncrease((PetscObject)vec));
1095: PetscFunctionReturn(PETSC_SUCCESS);
1096: }
1098: /*@
1099: VecLoad - Loads a vector that has been stored in binary or HDF5 format
1100: with `VecView()`.
1102: Collective
1104: Input Parameters:
1105: + vec - the newly loaded vector, this needs to have been created with `VecCreate()` or
1106: some related function before the call to `VecLoad()`.
1107: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()` or
1108: HDF5 file viewer, obtained from `PetscViewerHDF5Open()`
1110: Level: intermediate
1112: Notes:
1113: Defaults to the standard `VECSEQ` or `VECMPI`, if you want some other type of `Vec` call `VecSetFromOptions()`
1114: before calling this.
1116: The input file must contain the full global vector, as
1117: written by the routine `VecView()`.
1119: If the type or size of `vec` is not set before a call to `VecLoad()`, PETSc
1120: sets the type and the local and global sizes based on the vector it is reading in. If type and/or
1121: sizes are already set, then the same are used.
1123: If using the binary viewer and the blocksize of the vector is greater than one then you must provide a unique prefix to
1124: the vector with `PetscObjectSetOptionsPrefix`((`PetscObject`)vec,"uniqueprefix"); BEFORE calling `VecView()` on the
1125: vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
1126: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
1127: filename. If you copy the binary file, make sure you copy the associated .info file with it.
1129: If using HDF5, you must assign the `Vec` the same name as was used in the Vec
1130: that was stored in the file using `PetscObjectSetName(). Otherwise you will
1131: get the error message: "Cannot H5DOpen2() with `Vec` name NAMEOFOBJECT".
1133: If the HDF5 file contains a two dimensional array the first dimension is treated as the block size
1134: in loading the vector. Hence, for example, using MATLAB notation h5create('vector.dat','/Test_Vec',[27 1]);
1135: will load a vector of size 27 and block size 27 thus resulting in all 27 entries being on the first process of
1136: vectors communicator and the rest of the processes having zero entries
1138: Notes for advanced users when using the binary viewer:
1139: Most users should not need to know the details of the binary storage
1140: format, since `VecLoad()` and `VecView()` completely hide these details.
1141: But for anyone who's interested, the standard binary vector storage
1142: format is
1143: .vb
1144: PetscInt VEC_FILE_CLASSID
1145: PetscInt number of rows
1146: PetscScalar *values of all entries
1147: .ve
1149: In addition, PETSc automatically uses byte swapping to work on all machines; the files
1150: are written ALWAYS using big-endian ordering. On small-endian machines the numbers
1151: are converted to the small-endian format when they are read in from the file.
1152: See PetscBinaryRead() and PetscBinaryWrite() to see how this may be done.
1154: .seealso: [](ch_vectors), `Vec`, `PetscViewerBinaryOpen()`, `VecView()`, `MatLoad()`
1155: @*/
1156: PetscErrorCode VecLoad(Vec vec, PetscViewer viewer)
1157: {
1158: PetscBool isbinary, ishdf5, isadios, isexodusii, iscgns;
1159: PetscViewerFormat format;
1161: PetscFunctionBegin;
1164: PetscCheckSameComm(vec, 1, viewer, 2);
1165: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1166: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1167: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERCGNS, &iscgns));
1168: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERADIOS, &isadios));
1169: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWEREXODUSII, &isexodusii));
1170: PetscCheck(isbinary || ishdf5 || isadios || isexodusii || iscgns, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1172: PetscCall(VecSetErrorIfLocked(vec, 1));
1173: if (!((PetscObject)vec)->type_name && !vec->ops->create) PetscCall(VecSetType(vec, VECSTANDARD));
1174: PetscCall(PetscLogEventBegin(VEC_Load, viewer, 0, 0, 0));
1175: PetscCall(PetscViewerGetFormat(viewer, &format));
1176: if (format == PETSC_VIEWER_NATIVE && vec->ops->loadnative) {
1177: PetscUseTypeMethod(vec, loadnative, viewer);
1178: } else {
1179: PetscUseTypeMethod(vec, load, viewer);
1180: }
1181: PetscCall(PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0));
1182: PetscFunctionReturn(PETSC_SUCCESS);
1183: }
1185: /*@
1186: VecReciprocal - Replaces each component of a vector by its reciprocal.
1188: Logically Collective
1190: Input Parameter:
1191: . vec - the vector
1193: Output Parameter:
1194: . vec - the vector reciprocal
1196: Level: intermediate
1198: Note:
1199: Vector entries with value 0.0 are not changed
1201: .seealso: [](ch_vectors), `Vec`, `VecLog()`, `VecExp()`, `VecSqrtAbs()`
1202: @*/
1203: PetscErrorCode VecReciprocal(Vec vec)
1204: {
1205: PetscFunctionBegin;
1206: PetscCall(VecReciprocalAsync_Private(vec, NULL));
1207: PetscFunctionReturn(PETSC_SUCCESS);
1208: }
1210: /*@C
1211: VecSetOperation - Allows the user to override a particular vector operation.
1213: Logically Collective; No Fortran Support
1215: Input Parameters:
1216: + vec - The vector to modify
1217: . op - The name of the operation
1218: - f - The function that provides the operation.
1220: Level: advanced
1222: Example Usage:
1223: .vb
1224: // some new VecView() implementation, must have the same signature as the function it seeks
1225: // to replace
1226: PetscErrorCode UserVecView(Vec x, PetscViewer viewer)
1227: {
1228: PetscFunctionBeginUser;
1229: // ...
1230: PetscFunctionReturn(PETSC_SUCCESS);
1231: }
1233: // Create a VECMPI which has a pre-defined VecView() implementation
1234: VecCreateMPI(comm, n, N, &x);
1235: // Calls the VECMPI implementation for VecView()
1236: VecView(x, viewer);
1238: VecSetOperation(x, VECOP_VIEW, (PetscErrorCodeFn *)UserVecView);
1239: // Now calls UserVecView()
1240: VecView(x, viewer);
1241: .ve
1243: Notes:
1244: `f` may be `NULL` to remove the operation from `vec`. Depending on the operation this may be
1245: allowed, however some always expect a valid function. In these cases an error will be raised
1246: when calling the interface routine in question.
1248: See `VecOperation` for an up-to-date list of override-able operations. The operations listed
1249: there have the form `VECOP_<OPERATION>`, where `<OPERATION>` is the suffix (in all capital
1250: letters) of the public interface routine (e.g., `VecView()` -> `VECOP_VIEW`).
1252: Overriding a particular `Vec`'s operation has no affect on any other `Vec`s past, present,
1253: or future. The user should also note that overriding a method is "destructive"; the previous
1254: method is not retained in any way.
1256: Each function MUST return `PETSC_SUCCESS` on success and
1257: nonzero on failure.
1259: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecGetOperation()`, `MatSetOperation()`, `MatShellSetOperation()`
1260: @*/
1261: PetscErrorCode VecSetOperation(Vec vec, VecOperation op, PetscErrorCodeFn *f)
1262: {
1263: PetscFunctionBegin;
1265: if (op == VECOP_VIEW && !vec->ops->viewnative) {
1266: vec->ops->viewnative = vec->ops->view;
1267: } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
1268: vec->ops->loadnative = vec->ops->load;
1269: }
1270: ((PetscErrorCodeFn **)vec->ops)[(int)op] = f;
1271: PetscFunctionReturn(PETSC_SUCCESS);
1272: }
1274: /*@
1275: VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1276: used during the assembly process to store values that belong to
1277: other processors.
1279: Not Collective, different processes can have different size stashes
1281: Input Parameters:
1282: + vec - the vector
1283: . size - the initial size of the stash.
1284: - bsize - the initial size of the block-stash(if used).
1286: Options Database Keys:
1287: + -vecstash_initial_size <size> or <size0,size1,...sizep-1> - set initial size
1288: - -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1> - set initial block size
1290: Level: intermediate
1292: Notes:
1293: The block-stash is used for values set with `VecSetValuesBlocked()` while
1294: the stash is used for values set with `VecSetValues()`
1296: Run with the option -info and look for output of the form
1297: VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1298: to determine the appropriate value, MM, to use for size and
1299: VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1300: to determine the value, BMM to use for bsize
1302: PETSc attempts to smartly manage the stash size so there is little likelihood setting a
1303: a specific value here will affect performance
1305: .seealso: [](ch_vectors), `Vec`, `VecSetBlockSize()`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecStashView()`
1306: @*/
1307: PetscErrorCode VecStashSetInitialSize(Vec vec, PetscInt size, PetscInt bsize)
1308: {
1309: PetscFunctionBegin;
1311: PetscCall(VecStashSetInitialSize_Private(&vec->stash, size));
1312: PetscCall(VecStashSetInitialSize_Private(&vec->bstash, bsize));
1313: PetscFunctionReturn(PETSC_SUCCESS);
1314: }
1316: /*@
1317: VecSetRandom - Sets all components of a vector to random numbers.
1319: Logically Collective
1321: Input Parameters:
1322: + x - the vector
1323: - rctx - the random number context, formed by `PetscRandomCreate()`, or use `NULL` and it will create one internally.
1325: Output Parameter:
1326: . x - the vector
1328: Example of Usage:
1329: .vb
1330: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1331: VecSetRandom(x,rctx);
1332: PetscRandomDestroy(&rctx);
1333: .ve
1335: Level: intermediate
1337: .seealso: [](ch_vectors), `Vec`, `VecSet()`, `VecSetValues()`, `PetscRandomCreate()`, `PetscRandomDestroy()`
1338: @*/
1339: PetscErrorCode VecSetRandom(Vec x, PetscRandom rctx)
1340: {
1341: PetscRandom randObj = NULL;
1343: PetscFunctionBegin;
1347: VecCheckAssembled(x);
1348: PetscCall(VecSetErrorIfLocked(x, 1));
1350: if (!rctx) {
1351: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)x), &randObj));
1352: PetscCall(PetscRandomSetType(randObj, x->defaultrandtype));
1353: PetscCall(PetscRandomSetFromOptions(randObj));
1354: rctx = randObj;
1355: }
1357: PetscCall(PetscLogEventBegin(VEC_SetRandom, x, rctx, 0, 0));
1358: PetscUseTypeMethod(x, setrandom, rctx);
1359: PetscCall(PetscLogEventEnd(VEC_SetRandom, x, rctx, 0, 0));
1361: PetscCall(PetscRandomDestroy(&randObj));
1362: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1363: PetscFunctionReturn(PETSC_SUCCESS);
1364: }
1366: /*@
1367: VecZeroEntries - puts a `0.0` in each element of a vector
1369: Logically Collective
1371: Input Parameter:
1372: . vec - The vector
1374: Level: beginner
1376: Note:
1377: If the norm of the vector is known to be zero then this skips the unneeded zeroing process
1379: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecSetOptionsPrefix()`, `VecSet()`, `VecSetValues()`
1380: @*/
1381: PetscErrorCode VecZeroEntries(Vec vec)
1382: {
1383: PetscFunctionBegin;
1384: PetscCall(VecSet(vec, 0));
1385: PetscFunctionReturn(PETSC_SUCCESS);
1386: }
1388: /*
1389: VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
1390: processor and a PETSc MPI vector on more than one processor.
1392: Collective
1394: Input Parameter:
1395: . vec - The vector
1397: Level: intermediate
1399: .seealso: [](ch_vectors), `Vec`, `VecSetFromOptions()`, `VecSetType()`
1400: */
1401: static PetscErrorCode VecSetTypeFromOptions_Private(Vec vec, PetscOptionItems PetscOptionsObject)
1402: {
1403: PetscBool opt;
1404: VecType defaultType;
1405: char typeName[256];
1406: PetscMPIInt size;
1408: PetscFunctionBegin;
1409: if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1410: else {
1411: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
1412: if (size > 1) defaultType = VECMPI;
1413: else defaultType = VECSEQ;
1414: }
1416: PetscCall(VecRegisterAll());
1417: PetscCall(PetscOptionsFList("-vec_type", "Vector type", "VecSetType", VecList, defaultType, typeName, 256, &opt));
1418: if (opt) {
1419: PetscCall(VecSetType(vec, typeName));
1420: } else {
1421: PetscCall(VecSetType(vec, defaultType));
1422: }
1423: PetscFunctionReturn(PETSC_SUCCESS);
1424: }
1426: /*@
1427: VecSetFromOptions - Configures the vector from the options database.
1429: Collective
1431: Input Parameter:
1432: . vec - The vector
1434: Level: beginner
1436: Notes:
1437: To see all options, run your program with the -help option.
1439: Must be called after `VecCreate()` but before the vector is used.
1441: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecSetOptionsPrefix()`
1442: @*/
1443: PetscErrorCode VecSetFromOptions(Vec vec)
1444: {
1445: PetscBool flg;
1446: PetscInt bind_below = 0;
1448: PetscFunctionBegin;
1451: PetscObjectOptionsBegin((PetscObject)vec);
1452: /* Handle vector type options */
1453: PetscCall(VecSetTypeFromOptions_Private(vec, PetscOptionsObject));
1455: /* Handle specific vector options */
1456: PetscTryTypeMethod(vec, setfromoptions, PetscOptionsObject);
1458: /* Bind to CPU if below a user-specified size threshold.
1459: * This perhaps belongs in the options for the GPU Vec types, but VecBindToCPU() does nothing when called on non-GPU types,
1460: * and putting it here makes is more maintainable than duplicating this for all. */
1461: 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));
1462: if (flg && vec->map->n < bind_below) PetscCall(VecBindToCPU(vec, PETSC_TRUE));
1464: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1465: PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)vec, PetscOptionsObject));
1466: PetscOptionsEnd();
1467: PetscFunctionReturn(PETSC_SUCCESS);
1468: }
1470: /*@
1471: VecSetSizes - Sets the local and global sizes, and checks to determine compatibility of the sizes
1473: Collective
1475: Input Parameters:
1476: + v - the vector
1477: . n - the local size (or `PETSC_DECIDE` to have it set)
1478: - N - the global size (or `PETSC_DETERMINE` to have it set)
1480: Level: intermediate
1482: Notes:
1483: `N` cannot be `PETSC_DETERMINE` if `n` is `PETSC_DECIDE`
1485: If one processor calls this with `N` of `PETSC_DETERMINE` then all processors must, otherwise the program will hang.
1487: If `n` is not `PETSC_DECIDE`, then the value determines the `PetscLayout` of the vector and the ranges returned by
1488: `VecGetOwnershipRange()` and `VecGetOwnershipRanges()`
1490: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecCreateSeq()`, `VecCreateMPI()`, `VecGetSize()`, `PetscSplitOwnership()`, `PetscLayout`,
1491: `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`, `MatSetSizes()`
1492: @*/
1493: PetscErrorCode VecSetSizes(Vec v, PetscInt n, PetscInt N)
1494: {
1495: PetscFunctionBegin;
1497: if (N >= 0) {
1499: PetscCheck(n <= N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local size %" PetscInt_FMT " cannot be larger than global size %" PetscInt_FMT, n, N);
1500: }
1501: 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,
1502: v->map->n, v->map->N);
1503: v->map->n = n;
1504: v->map->N = N;
1505: PetscTryTypeMethod(v, create);
1506: v->ops->create = NULL;
1507: PetscFunctionReturn(PETSC_SUCCESS);
1508: }
1510: /*@
1511: VecSetBlockSize - Sets the block size for future calls to `VecSetValuesBlocked()`
1512: and `VecSetValuesBlockedLocal()`.
1514: Logically Collective
1516: Input Parameters:
1517: + v - the vector
1518: - bs - the blocksize
1520: Level: advanced
1522: Note:
1523: All vectors obtained by `VecDuplicate()` inherit the same blocksize.
1525: Vectors obtained with `DMCreateGlobalVector()` and `DMCreateLocalVector()` generally already have a blocksize set based on the state of the `DM`
1527: .seealso: [](ch_vectors), `Vec`, `VecSetValuesBlocked()`, `VecSetLocalToGlobalMapping()`, `VecGetBlockSize()`
1528: @*/
1529: PetscErrorCode VecSetBlockSize(Vec v, PetscInt bs)
1530: {
1531: PetscFunctionBegin;
1534: PetscCall(PetscLayoutSetBlockSize(v->map, bs));
1535: v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1536: PetscFunctionReturn(PETSC_SUCCESS);
1537: }
1539: /*@
1540: VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for `VecSetValuesBlocked()`
1541: and `VecSetValuesBlockedLocal()`.
1543: Not Collective
1545: Input Parameter:
1546: . v - the vector
1548: Output Parameter:
1549: . bs - the blocksize
1551: Level: advanced
1553: Note:
1554: All vectors obtained by `VecDuplicate()` inherit the same blocksize.
1556: .seealso: [](ch_vectors), `Vec`, `VecSetValuesBlocked()`, `VecSetLocalToGlobalMapping()`, `VecSetBlockSize()`
1557: @*/
1558: PetscErrorCode VecGetBlockSize(Vec v, PetscInt *bs)
1559: {
1560: PetscFunctionBegin;
1562: PetscAssertPointer(bs, 2);
1563: PetscCall(PetscLayoutGetBlockSize(v->map, bs));
1564: PetscFunctionReturn(PETSC_SUCCESS);
1565: }
1567: /*@
1568: VecSetOptionsPrefix - Sets the prefix used for searching for all
1569: `Vec` options in the database.
1571: Logically Collective
1573: Input Parameters:
1574: + v - the `Vec` context
1575: - prefix - the prefix to prepend to all option names
1577: Level: advanced
1579: Note:
1580: A hyphen (-) must NOT be given at the beginning of the prefix name.
1581: The first character of all runtime options is AUTOMATICALLY the hyphen.
1583: .seealso: [](ch_vectors), `Vec`, `VecSetFromOptions()`
1584: @*/
1585: PetscErrorCode VecSetOptionsPrefix(Vec v, const char prefix[])
1586: {
1587: PetscFunctionBegin;
1589: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)v, prefix));
1590: PetscFunctionReturn(PETSC_SUCCESS);
1591: }
1593: /*@
1594: VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1595: `Vec` options in the database.
1597: Logically Collective
1599: Input Parameters:
1600: + v - the `Vec` context
1601: - prefix - the prefix to prepend to all option names
1603: Level: advanced
1605: Note:
1606: A hyphen (-) must NOT be given at the beginning of the prefix name.
1607: The first character of all runtime options is AUTOMATICALLY the hyphen.
1609: .seealso: [](ch_vectors), `Vec`, `VecGetOptionsPrefix()`
1610: @*/
1611: PetscErrorCode VecAppendOptionsPrefix(Vec v, const char prefix[])
1612: {
1613: PetscFunctionBegin;
1615: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)v, prefix));
1616: PetscFunctionReturn(PETSC_SUCCESS);
1617: }
1619: /*@
1620: VecGetOptionsPrefix - Sets the prefix used for searching for all
1621: Vec options in the database.
1623: Not Collective
1625: Input Parameter:
1626: . v - the `Vec` context
1628: Output Parameter:
1629: . prefix - pointer to the prefix string used
1631: Level: advanced
1633: .seealso: [](ch_vectors), `Vec`, `VecAppendOptionsPrefix()`
1634: @*/
1635: PetscErrorCode VecGetOptionsPrefix(Vec v, const char *prefix[])
1636: {
1637: PetscFunctionBegin;
1639: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)v, prefix));
1640: PetscFunctionReturn(PETSC_SUCCESS);
1641: }
1643: /*@C
1644: VecGetState - Gets the state of a `Vec`.
1646: Not Collective
1648: Input Parameter:
1649: . v - the `Vec` context
1651: Output Parameter:
1652: . state - the object state
1654: Level: advanced
1656: Note:
1657: Object state is an integer which gets increased every time
1658: the object is changed. By saving and later querying the object state
1659: one can determine whether information about the object is still current.
1661: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `PetscObjectStateGet()`
1662: @*/
1663: PetscErrorCode VecGetState(Vec v, PetscObjectState *state)
1664: {
1665: PetscFunctionBegin;
1667: PetscAssertPointer(state, 2);
1668: PetscCall(PetscObjectStateGet((PetscObject)v, state));
1669: PetscFunctionReturn(PETSC_SUCCESS);
1670: }
1672: /*@
1673: VecSetUp - Sets up the internal vector data structures for the later use.
1675: Collective
1677: Input Parameter:
1678: . v - the `Vec` context
1680: Level: advanced
1682: Notes:
1683: For basic use of the `Vec` classes the user need not explicitly call
1684: `VecSetUp()`, since these actions will happen automatically.
1686: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecDestroy()`
1687: @*/
1688: PetscErrorCode VecSetUp(Vec v)
1689: {
1690: PetscMPIInt size;
1692: PetscFunctionBegin;
1694: PetscCheck(v->map->n >= 0 || v->map->N >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Sizes not set");
1695: if (!((PetscObject)v)->type_name) {
1696: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1697: if (size == 1) {
1698: PetscCall(VecSetType(v, VECSEQ));
1699: } else {
1700: PetscCall(VecSetType(v, VECMPI));
1701: }
1702: }
1703: PetscFunctionReturn(PETSC_SUCCESS);
1704: }
1706: /*
1707: These currently expose the PetscScalar/PetscReal in updating the
1708: cached norm. If we push those down into the implementation these
1709: will become independent of PetscScalar/PetscReal
1710: */
1712: PetscErrorCode VecCopyAsync_Private(Vec x, Vec y, PetscDeviceContext dctx)
1713: {
1714: PetscBool flgs[4];
1715: PetscReal norms[4] = {0.0, 0.0, 0.0, 0.0};
1717: PetscFunctionBegin;
1722: if (x == y) PetscFunctionReturn(PETSC_SUCCESS);
1723: VecCheckSameLocalSize(x, 1, y, 2);
1724: VecCheckAssembled(x);
1725: PetscCall(VecSetErrorIfLocked(y, 2));
1727: #if !defined(PETSC_USE_MIXED_PRECISION)
1728: for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
1729: #endif
1731: PetscCall(PetscLogEventBegin(VEC_Copy, x, y, 0, 0));
1732: #if defined(PETSC_USE_MIXED_PRECISION)
1733: extern PetscErrorCode VecGetArray(Vec, double **);
1734: extern PetscErrorCode VecRestoreArray(Vec, double **);
1735: extern PetscErrorCode VecGetArray(Vec, float **);
1736: extern PetscErrorCode VecRestoreArray(Vec, float **);
1737: extern PetscErrorCode VecGetArrayRead(Vec, const double **);
1738: extern PetscErrorCode VecRestoreArrayRead(Vec, const double **);
1739: extern PetscErrorCode VecGetArrayRead(Vec, const float **);
1740: extern PetscErrorCode VecRestoreArrayRead(Vec, const float **);
1741: if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1742: PetscInt i, n;
1743: const float *xx;
1744: double *yy;
1745: PetscCall(VecGetArrayRead(x, &xx));
1746: PetscCall(VecGetArray(y, &yy));
1747: PetscCall(VecGetLocalSize(x, &n));
1748: for (i = 0; i < n; i++) yy[i] = xx[i];
1749: PetscCall(VecRestoreArrayRead(x, &xx));
1750: PetscCall(VecRestoreArray(y, &yy));
1751: } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1752: PetscInt i, n;
1753: float *yy;
1754: const double *xx;
1755: PetscCall(VecGetArrayRead(x, &xx));
1756: PetscCall(VecGetArray(y, &yy));
1757: PetscCall(VecGetLocalSize(x, &n));
1758: for (i = 0; i < n; i++) yy[i] = (float)xx[i];
1759: PetscCall(VecRestoreArrayRead(x, &xx));
1760: PetscCall(VecRestoreArray(y, &yy));
1761: } else PetscUseTypeMethod(x, copy, y);
1762: #else
1763: VecMethodDispatch(x, dctx, VecAsyncFnName(Copy), copy, (Vec, Vec, PetscDeviceContext), y);
1764: #endif
1766: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1767: #if !defined(PETSC_USE_MIXED_PRECISION)
1768: for (PetscInt i = 0; i < 4; i++) {
1769: if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)y, NormIds[i], norms[i]));
1770: }
1771: #endif
1773: PetscCall(PetscLogEventEnd(VEC_Copy, x, y, 0, 0));
1774: PetscFunctionReturn(PETSC_SUCCESS);
1775: }
1777: /*@
1778: VecCopy - Copies a vector `y = x`
1780: Logically Collective
1782: Input Parameter:
1783: . x - the vector
1785: Output Parameter:
1786: . y - the copy
1788: Level: beginner
1790: Note:
1791: For default parallel PETSc vectors, both `x` and `y` must be distributed in
1792: the same manner; local copies are done.
1794: Developer Notes:
1795: `PetscCheckSameTypeAndComm`(x,1,y,2) is not used on these vectors because we allow one
1796: of the vectors to be sequential and one to be parallel so long as both have the same
1797: local sizes. This is used in some internal functions in PETSc.
1799: .seealso: [](ch_vectors), `Vec`, `VecDuplicate()`
1800: @*/
1801: PetscErrorCode VecCopy(Vec x, Vec y)
1802: {
1803: PetscFunctionBegin;
1804: PetscCall(VecCopyAsync_Private(x, y, NULL));
1805: PetscFunctionReturn(PETSC_SUCCESS);
1806: }
1808: PetscErrorCode VecSwapAsync_Private(Vec x, Vec y, PetscDeviceContext dctx)
1809: {
1810: PetscReal normxs[4], normys[4];
1811: PetscBool flgxs[4], flgys[4];
1813: PetscFunctionBegin;
1818: PetscCheckSameTypeAndComm(x, 1, y, 2);
1819: VecCheckSameSize(x, 1, y, 2);
1820: VecCheckAssembled(x);
1821: VecCheckAssembled(y);
1822: PetscCall(VecSetErrorIfLocked(x, 1));
1823: PetscCall(VecSetErrorIfLocked(y, 2));
1825: for (PetscInt i = 0; i < 4; i++) {
1826: PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], normxs[i], flgxs[i]));
1827: PetscCall(PetscObjectComposedDataGetReal((PetscObject)y, NormIds[i], normys[i], flgys[i]));
1828: }
1830: PetscCall(PetscLogEventBegin(VEC_Swap, x, y, 0, 0));
1831: VecMethodDispatch(x, dctx, VecAsyncFnName(Swap), swap, (Vec, Vec, PetscDeviceContext), y);
1832: PetscCall(PetscLogEventEnd(VEC_Swap, x, y, 0, 0));
1834: PetscCall(PetscObjectStateIncrease((PetscObject)x));
1835: PetscCall(PetscObjectStateIncrease((PetscObject)y));
1836: for (PetscInt i = 0; i < 4; i++) {
1837: if (flgxs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)y, NormIds[i], normxs[i]));
1838: if (flgys[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], normys[i]));
1839: }
1840: PetscFunctionReturn(PETSC_SUCCESS);
1841: }
1842: /*@
1843: VecSwap - Swaps the values between two vectors, `x` and `y`.
1845: Logically Collective
1847: Input Parameters:
1848: + x - the first vector
1849: - y - the second vector
1851: Level: advanced
1853: .seealso: [](ch_vectors), `Vec`, `VecSet()`
1854: @*/
1855: PetscErrorCode VecSwap(Vec x, Vec y)
1856: {
1857: PetscFunctionBegin;
1858: PetscCall(VecSwapAsync_Private(x, y, NULL));
1859: PetscFunctionReturn(PETSC_SUCCESS);
1860: }
1862: /*@
1863: VecStashViewFromOptions - Processes command line options to determine if/how a `VecStash` object is to be viewed.
1865: Collective
1867: Input Parameters:
1868: + obj - the `Vec` containing a stash
1869: . bobj - optional other object that provides the prefix
1870: - optionname - option to activate viewing
1872: Level: intermediate
1874: Developer Notes:
1875: This cannot use `PetscObjectViewFromOptions()` because it takes a `Vec` as an argument but does not use `VecView()`
1877: .seealso: [](ch_vectors), `Vec`, `VecStashSetInitialSize()`
1878: @*/
1879: PetscErrorCode VecStashViewFromOptions(Vec obj, PetscObject bobj, const char optionname[])
1880: {
1881: PetscViewer viewer;
1882: PetscBool flg;
1883: PetscViewerFormat format;
1884: char *prefix;
1886: PetscFunctionBegin;
1887: prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
1888: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)obj), ((PetscObject)obj)->options, prefix, optionname, &viewer, &format, &flg));
1889: if (flg) {
1890: PetscCall(PetscViewerPushFormat(viewer, format));
1891: PetscCall(VecStashView(obj, viewer));
1892: PetscCall(PetscViewerPopFormat(viewer));
1893: PetscCall(PetscViewerDestroy(&viewer));
1894: }
1895: PetscFunctionReturn(PETSC_SUCCESS);
1896: }
1898: /*@
1899: VecStashView - Prints the entries in the vector stash and block stash.
1901: Collective
1903: Input Parameters:
1904: + v - the vector
1905: - viewer - the viewer
1907: Level: advanced
1909: .seealso: [](ch_vectors), `Vec`, `VecSetBlockSize()`, `VecSetValues()`, `VecSetValuesBlocked()`
1910: @*/
1911: PetscErrorCode VecStashView(Vec v, PetscViewer viewer)
1912: {
1913: PetscMPIInt rank;
1914: PetscInt i, j;
1915: PetscBool match;
1916: VecStash *s;
1917: PetscScalar val;
1919: PetscFunctionBegin;
1922: PetscCheckSameComm(v, 1, viewer, 2);
1924: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &match));
1925: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_SUP, "Stash viewer only works with ASCII viewer not %s", ((PetscObject)v)->type_name);
1926: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1927: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank));
1928: s = &v->bstash;
1930: /* print block stash */
1931: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
1932: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Vector Block stash size %" PetscInt_FMT " block size %" PetscInt_FMT "\n", rank, s->n, s->bs));
1933: for (i = 0; i < s->n; i++) {
1934: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " ", rank, s->idx[i]));
1935: for (j = 0; j < s->bs; j++) {
1936: val = s->array[i * s->bs + j];
1937: #if defined(PETSC_USE_COMPLEX)
1938: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "(%18.16e %18.16e) ", (double)PetscRealPart(val), (double)PetscImaginaryPart(val)));
1939: #else
1940: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%18.16e ", (double)val));
1941: #endif
1942: }
1943: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
1944: }
1945: PetscCall(PetscViewerFlush(viewer));
1947: s = &v->stash;
1949: /* print basic stash */
1950: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Vector stash size %" PetscInt_FMT "\n", rank, s->n));
1951: for (i = 0; i < s->n; i++) {
1952: val = s->array[i];
1953: #if defined(PETSC_USE_COMPLEX)
1954: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " (%18.16e %18.16e) ", rank, s->idx[i], (double)PetscRealPart(val), (double)PetscImaginaryPart(val)));
1955: #else
1956: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " %18.16e\n", rank, s->idx[i], (double)val));
1957: #endif
1958: }
1959: PetscCall(PetscViewerFlush(viewer));
1960: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1961: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1962: PetscFunctionReturn(PETSC_SUCCESS);
1963: }
1965: PetscErrorCode PetscOptionsGetVec(PetscOptions options, const char prefix[], const char key[], Vec v, PetscBool *set)
1966: {
1967: PetscInt i, N, rstart, rend;
1968: PetscScalar *xx;
1969: PetscReal *xreal;
1970: PetscBool iset;
1972: PetscFunctionBegin;
1973: PetscCall(VecGetOwnershipRange(v, &rstart, &rend));
1974: PetscCall(VecGetSize(v, &N));
1975: PetscCall(PetscCalloc1(N, &xreal));
1976: PetscCall(PetscOptionsGetRealArray(options, prefix, key, xreal, &N, &iset));
1977: if (iset) {
1978: PetscCall(VecGetArray(v, &xx));
1979: for (i = rstart; i < rend; i++) xx[i - rstart] = xreal[i];
1980: PetscCall(VecRestoreArray(v, &xx));
1981: }
1982: PetscCall(PetscFree(xreal));
1983: if (set) *set = iset;
1984: PetscFunctionReturn(PETSC_SUCCESS);
1985: }
1987: /*@
1988: VecGetLayout - get `PetscLayout` describing a vector layout
1990: Not Collective
1992: Input Parameter:
1993: . x - the vector
1995: Output Parameter:
1996: . map - the layout
1998: Level: developer
2000: Note:
2001: The layout determines what vector elements are contained on each MPI process
2003: .seealso: [](ch_vectors), `PetscLayout`, `Vec`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
2004: @*/
2005: PetscErrorCode VecGetLayout(Vec x, PetscLayout *map)
2006: {
2007: PetscFunctionBegin;
2009: PetscAssertPointer(map, 2);
2010: *map = x->map;
2011: PetscFunctionReturn(PETSC_SUCCESS);
2012: }
2014: /*@
2015: VecSetLayout - set `PetscLayout` describing vector layout
2017: Not Collective
2019: Input Parameters:
2020: + x - the vector
2021: - map - the layout
2023: Level: developer
2025: Note:
2026: It is normally only valid to replace the layout with a layout known to be equivalent.
2028: .seealso: [](ch_vectors), `Vec`, `PetscLayout`, `VecGetLayout()`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
2029: @*/
2030: PetscErrorCode VecSetLayout(Vec x, PetscLayout map)
2031: {
2032: PetscFunctionBegin;
2034: PetscCall(PetscLayoutReference(map, &x->map));
2035: PetscFunctionReturn(PETSC_SUCCESS);
2036: }
2038: /*@
2039: VecFlag - set infinity into the local part of the vector on any subset of MPI processes
2041: Logically Collective
2043: Input Parameters:
2044: + xin - the vector, can be `NULL` but only if on all processes
2045: - flg - indicates if this processes portion of the vector should be set to infinity
2047: Level: developer
2049: Note:
2050: This removes the values from the vector norm cache for all processes by calling `PetscObjectIncrease()`.
2052: This is used for any subset of MPI processes to indicate an failure in a solver, after the next use of `VecNorm()` if
2053: `KSPCheckNorm()` detects an infinity and at least one of the MPI processes has a not converged reason then the `KSP`
2054: object collectively is labeled as not converged.
2056: .seealso: [](ch_vectors), `Vec`, `PetscLayout`, `VecGetLayout()`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
2057: @*/
2058: PetscErrorCode VecFlag(Vec xin, PetscInt flg)
2059: {
2060: // MSVC gives "divide by zero" error at compile time - so declare as volatile to skip this check.
2061: volatile PetscReal one = 1.0, zero = 0.0;
2062: PetscScalar inf;
2064: PetscFunctionBegin;
2065: if (!xin) PetscFunctionReturn(PETSC_SUCCESS);
2067: PetscCall(PetscObjectStateIncrease((PetscObject)xin));
2068: if (flg) {
2069: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2070: inf = one / zero;
2071: PetscCall(PetscFPTrapPop());
2072: if (xin->ops->set) {
2073: PetscUseTypeMethod(xin, set, inf);
2074: } else {
2075: PetscInt n;
2076: PetscScalar *xx;
2078: PetscCall(VecGetLocalSize(xin, &n));
2079: PetscCall(VecGetArrayWrite(xin, &xx));
2080: for (PetscInt i = 0; i < n; ++i) xx[i] = inf;
2081: PetscCall(VecRestoreArrayWrite(xin, &xx));
2082: }
2083: }
2084: PetscFunctionReturn(PETSC_SUCCESS);
2085: }
2087: /*@
2088: VecSetInf - set infinity into the local part of the vector
2090: Not Collective
2092: Input Parameters:
2093: . xin - the vector
2095: Level: developer
2097: Note:
2098: Deprecated, see `VecFlag()`
2099: This is used for any subset of MPI processes to indicate an failure in a solver, after the next use of `VecNorm()` if
2100: `KSPCheckNorm()` detects an infinity and at least one of the MPI processes has a not converged reason then the `KSP`
2101: object collectively is labeled as not converged.
2103: This cannot be called if `xin` has a cached norm available
2105: .seealso: [](ch_vectors), `VecFlag()`, `Vec`, `PetscLayout`, `VecGetLayout()`, `VecGetSize()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
2106: @*/
2107: PetscErrorCode VecSetInf(Vec xin)
2108: {
2109: // MSVC gives "divide by zero" error at compile time - so declare as volatile to skip this check.
2110: volatile PetscReal one = 1.0, zero = 0.0;
2111: PetscScalar inf;
2112: PetscBool flg;
2114: PetscFunctionBegin;
2115: PetscCall(VecNormAvailable(xin, NORM_2, &flg, NULL));
2116: PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot call VecSetInf() if the vector has a cached norm");
2117: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2118: inf = one / zero;
2119: PetscCall(PetscFPTrapPop());
2120: if (xin->ops->set) {
2121: PetscUseTypeMethod(xin, set, inf);
2122: } else {
2123: PetscInt n;
2124: PetscScalar *xx;
2126: PetscCall(VecGetLocalSize(xin, &n));
2127: PetscCall(VecGetArrayWrite(xin, &xx));
2128: for (PetscInt i = 0; i < n; ++i) xx[i] = inf;
2129: PetscCall(VecRestoreArrayWrite(xin, &xx));
2130: }
2131: PetscFunctionReturn(PETSC_SUCCESS);
2132: }
2134: /*@
2135: VecBindToCPU - marks a vector to temporarily stay on the CPU and perform computations on the CPU
2137: Logically collective
2139: Input Parameters:
2140: + v - the vector
2141: - flg - bind to the CPU if value of `PETSC_TRUE`
2143: Level: intermediate
2145: .seealso: [](ch_vectors), `Vec`, `VecBoundToCPU()`
2146: @*/
2147: PetscErrorCode VecBindToCPU(Vec v, PetscBool flg)
2148: {
2149: PetscFunctionBegin;
2152: #if defined(PETSC_HAVE_DEVICE)
2153: if (v->boundtocpu == flg) PetscFunctionReturn(PETSC_SUCCESS);
2154: v->boundtocpu = flg;
2155: PetscTryTypeMethod(v, bindtocpu, flg);
2156: #endif
2157: PetscFunctionReturn(PETSC_SUCCESS);
2158: }
2160: /*@
2161: VecBoundToCPU - query if a vector is bound to the CPU
2163: Not collective
2165: Input Parameter:
2166: . v - the vector
2168: Output Parameter:
2169: . flg - the logical flag
2171: Level: intermediate
2173: .seealso: [](ch_vectors), `Vec`, `VecBindToCPU()`
2174: @*/
2175: PetscErrorCode VecBoundToCPU(Vec v, PetscBool *flg)
2176: {
2177: PetscFunctionBegin;
2179: PetscAssertPointer(flg, 2);
2180: #if defined(PETSC_HAVE_DEVICE)
2181: *flg = v->boundtocpu;
2182: #else
2183: *flg = PETSC_TRUE;
2184: #endif
2185: PetscFunctionReturn(PETSC_SUCCESS);
2186: }
2188: /*@
2189: VecSetBindingPropagates - Sets whether the state of being bound to the CPU for a GPU vector type propagates to child and some other associated objects
2191: Input Parameters:
2192: + v - the vector
2193: - flg - flag indicating whether the boundtocpu flag should be propagated
2195: Level: developer
2197: Notes:
2198: 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.
2199: The created vectors will also have their bindingpropagates flag set to true.
2201: Developer Notes:
2202: If a `DMDA` has the `-dm_bind_below option` set to true, then vectors created by `DMCreateGlobalVector()` will have `VecSetBindingPropagates()` called on them to
2203: set their bindingpropagates flag to true.
2205: .seealso: [](ch_vectors), `Vec`, `MatSetBindingPropagates()`, `VecGetBindingPropagates()`
2206: @*/
2207: PetscErrorCode VecSetBindingPropagates(Vec v, PetscBool flg)
2208: {
2209: PetscFunctionBegin;
2211: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2212: v->bindingpropagates = flg;
2213: #endif
2214: PetscFunctionReturn(PETSC_SUCCESS);
2215: }
2217: /*@
2218: VecGetBindingPropagates - Gets whether the state of being bound to the CPU for a GPU vector type propagates to child and some other associated objects
2220: Input Parameter:
2221: . v - the vector
2223: Output Parameter:
2224: . flg - flag indicating whether the boundtocpu flag will be propagated
2226: Level: developer
2228: .seealso: [](ch_vectors), `Vec`, `VecSetBindingPropagates()`
2229: @*/
2230: PetscErrorCode VecGetBindingPropagates(Vec v, PetscBool *flg)
2231: {
2232: PetscFunctionBegin;
2234: PetscAssertPointer(flg, 2);
2235: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2236: *flg = v->bindingpropagates;
2237: #else
2238: *flg = PETSC_FALSE;
2239: #endif
2240: PetscFunctionReturn(PETSC_SUCCESS);
2241: }
2243: /*@C
2244: VecSetPinnedMemoryMin - Set the minimum data size for which pinned memory will be used for host (CPU) allocations.
2246: Logically Collective
2248: Input Parameters:
2249: + v - the vector
2250: - mbytes - minimum data size in bytes
2252: Options Database Key:
2253: . -vec_pinned_memory_min <size> - minimum size (in bytes) for an allocation to use pinned memory on host.
2255: Level: developer
2257: Note:
2258: Specifying -1 ensures that pinned memory will never be used.
2260: .seealso: [](ch_vectors), `Vec`, `VecGetPinnedMemoryMin()`
2261: @*/
2262: PetscErrorCode VecSetPinnedMemoryMin(Vec v, size_t mbytes)
2263: {
2264: PetscFunctionBegin;
2266: #if PetscDefined(HAVE_DEVICE)
2267: v->minimum_bytes_pinned_memory = mbytes;
2268: #endif
2269: PetscFunctionReturn(PETSC_SUCCESS);
2270: }
2272: /*@C
2273: VecGetPinnedMemoryMin - Get the minimum data size for which pinned memory will be used for host (CPU) allocations.
2275: Logically Collective
2277: Input Parameter:
2278: . v - the vector
2280: Output Parameter:
2281: . mbytes - minimum data size in bytes
2283: Level: developer
2285: .seealso: [](ch_vectors), `Vec`, `VecSetPinnedMemoryMin()`
2286: @*/
2287: PetscErrorCode VecGetPinnedMemoryMin(Vec v, size_t *mbytes)
2288: {
2289: PetscFunctionBegin;
2291: PetscAssertPointer(mbytes, 2);
2292: #if PetscDefined(HAVE_DEVICE)
2293: *mbytes = v->minimum_bytes_pinned_memory;
2294: #endif
2295: PetscFunctionReturn(PETSC_SUCCESS);
2296: }
2298: /*@
2299: VecGetOffloadMask - Get the offload mask of a `Vec`
2301: Not Collective
2303: Input Parameter:
2304: . v - the vector
2306: Output Parameter:
2307: . mask - corresponding `PetscOffloadMask` enum value.
2309: Level: intermediate
2311: .seealso: [](ch_vectors), `Vec`, `VecCreateSeqCUDA()`, `VecCreateSeqViennaCL()`, `VecGetArray()`, `VecGetType()`
2312: @*/
2313: PetscErrorCode VecGetOffloadMask(Vec v, PetscOffloadMask *mask)
2314: {
2315: PetscFunctionBegin;
2317: PetscAssertPointer(mask, 2);
2318: *mask = v->offloadmask;
2319: PetscFunctionReturn(PETSC_SUCCESS);
2320: }
2322: #if !defined(PETSC_HAVE_VIENNACL)
2323: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLContext(Vec v, PETSC_UINTPTR_T *ctx)
2324: {
2325: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_context");
2326: }
2328: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLQueue(Vec v, PETSC_UINTPTR_T *queue)
2329: {
2330: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_command_queue");
2331: }
2333: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMem(Vec v, PETSC_UINTPTR_T *queue)
2334: {
2335: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2336: }
2338: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMemRead(Vec v, PETSC_UINTPTR_T *queue)
2339: {
2340: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2341: }
2343: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMemWrite(Vec v, PETSC_UINTPTR_T *queue)
2344: {
2345: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2346: }
2348: PETSC_EXTERN PetscErrorCode VecViennaCLRestoreCLMemWrite(Vec v)
2349: {
2350: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to restore a Vec's cl_mem");
2351: }
2352: #endif
2354: 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)
2355: {
2356: const PetscScalar *u, *y;
2357: const PetscScalar *atola = NULL, *rtola = NULL, *erra = NULL;
2358: PetscInt n, n_loc = 0, na_loc = 0, nr_loc = 0;
2359: PetscReal nrm = 0, nrma = 0, nrmr = 0, err_loc[6];
2361: PetscFunctionBegin;
2362: #define SkipSmallValue(a, b, tol) \
2363: if (PetscAbsScalar(a) < tol || PetscAbsScalar(b) < tol) continue
2365: PetscCall(VecGetLocalSize(U, &n));
2366: PetscCall(VecGetArrayRead(U, &u));
2367: PetscCall(VecGetArrayRead(Y, &y));
2368: if (E) PetscCall(VecGetArrayRead(E, &erra));
2369: if (vatol) PetscCall(VecGetArrayRead(vatol, &atola));
2370: if (vrtol) PetscCall(VecGetArrayRead(vrtol, &rtola));
2371: for (PetscInt i = 0; i < n; i++) {
2372: PetscReal err, tol, tola, tolr;
2374: SkipSmallValue(y[i], u[i], ignore_max);
2375: atol = atola ? PetscRealPart(atola[i]) : atol;
2376: rtol = rtola ? PetscRealPart(rtola[i]) : rtol;
2377: err = erra ? PetscAbsScalar(erra[i]) : PetscAbsScalar(y[i] - u[i]);
2378: tola = atol;
2379: tolr = rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
2380: tol = tola + tolr;
2381: if (tola > 0.) {
2382: if (wnormtype == NORM_INFINITY) nrma = PetscMax(nrma, err / tola);
2383: else nrma += PetscSqr(err / tola);
2384: na_loc++;
2385: }
2386: if (tolr > 0.) {
2387: if (wnormtype == NORM_INFINITY) nrmr = PetscMax(nrmr, err / tolr);
2388: else nrmr += PetscSqr(err / tolr);
2389: nr_loc++;
2390: }
2391: if (tol > 0.) {
2392: if (wnormtype == NORM_INFINITY) nrm = PetscMax(nrm, err / tol);
2393: else nrm += PetscSqr(err / tol);
2394: n_loc++;
2395: }
2396: }
2397: if (E) PetscCall(VecRestoreArrayRead(E, &erra));
2398: if (vatol) PetscCall(VecRestoreArrayRead(vatol, &atola));
2399: if (vrtol) PetscCall(VecRestoreArrayRead(vrtol, &rtola));
2400: PetscCall(VecRestoreArrayRead(U, &u));
2401: PetscCall(VecRestoreArrayRead(Y, &y));
2402: #undef SkipSmallValue
2404: err_loc[0] = nrm;
2405: err_loc[1] = nrma;
2406: err_loc[2] = nrmr;
2407: err_loc[3] = (PetscReal)n_loc;
2408: err_loc[4] = (PetscReal)na_loc;
2409: err_loc[5] = (PetscReal)nr_loc;
2410: if (wnormtype == NORM_2) {
2411: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, err_loc, 6, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)U)));
2412: } else {
2413: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, err_loc, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)U)));
2414: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, err_loc + 3, 3, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)U)));
2415: }
2416: if (wnormtype == NORM_2) {
2417: *norm = PetscSqrtReal(err_loc[0]);
2418: *norma = PetscSqrtReal(err_loc[1]);
2419: *normr = PetscSqrtReal(err_loc[2]);
2420: } else {
2421: *norm = err_loc[0];
2422: *norma = err_loc[1];
2423: *normr = err_loc[2];
2424: }
2425: *norm_loc = (PetscInt)err_loc[3];
2426: *norma_loc = (PetscInt)err_loc[4];
2427: *normr_loc = (PetscInt)err_loc[5];
2428: PetscFunctionReturn(PETSC_SUCCESS);
2429: }
2431: /*@
2432: VecErrorWeightedNorms - compute a weighted norm of the difference between two vectors
2434: Collective
2436: Input Parameters:
2437: + U - first vector to be compared
2438: . Y - second vector to be compared
2439: . E - optional third vector representing the error (if not provided, the error is ||U-Y||)
2440: . wnormtype - norm type
2441: . atol - scalar for absolute tolerance
2442: . vatol - vector representing per-entry absolute tolerances (can be ``NULL``)
2443: . rtol - scalar for relative tolerance
2444: . vrtol - vector representing per-entry relative tolerances (can be ``NULL``)
2445: - ignore_max - ignore values smaller than this value in absolute terms.
2447: Output Parameters:
2448: + norm - weighted norm
2449: . norm_loc - number of vector locations used for the weighted norm
2450: . norma - weighted norm based on the absolute tolerance
2451: . norma_loc - number of vector locations used for the absolute weighted norm
2452: . normr - weighted norm based on the relative tolerance
2453: - normr_loc - number of vector locations used for the relative weighted norm
2455: Level: developer
2457: Notes:
2458: This is primarily used for computing weighted local truncation errors in ``TS``.
2460: .seealso: [](ch_vectors), `Vec`, `NormType`, `TSErrorWeightedNorm()`, `TSErrorWeightedENorm()`
2461: @*/
2462: 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)
2463: {
2464: PetscFunctionBegin;
2469: if (E) {
2472: }
2475: if (vatol) {
2478: }
2480: if (vrtol) {
2483: }
2485: PetscAssertPointer(norm, 10);
2486: PetscAssertPointer(norm_loc, 11);
2487: PetscAssertPointer(norma, 12);
2488: PetscAssertPointer(norma_loc, 13);
2489: PetscAssertPointer(normr, 14);
2490: PetscAssertPointer(normr_loc, 15);
2491: PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)U), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);
2493: /* There are potentially 5 vectors involved, some of them may happen to be of different type or bound to cpu.
2494: Here we check that they all implement the same operation and call it if so.
2495: Otherwise, we call the _Basic implementation that always works (provided VecGetArrayRead is implemented). */
2496: PetscBool sameop = (PetscBool)(U->ops->errorwnorm && U->ops->errorwnorm == Y->ops->errorwnorm);
2497: if (sameop && E) sameop = (PetscBool)(U->ops->errorwnorm == E->ops->errorwnorm);
2498: if (sameop && vatol) sameop = (PetscBool)(U->ops->errorwnorm == vatol->ops->errorwnorm);
2499: if (sameop && vrtol) sameop = (PetscBool)(U->ops->errorwnorm == vrtol->ops->errorwnorm);
2500: if (sameop) PetscUseTypeMethod(U, errorwnorm, Y, E, wnormtype, atol, vatol, rtol, vrtol, ignore_max, norm, norm_loc, norma, norma_loc, normr, normr_loc);
2501: else PetscCall(VecErrorWeightedNorms_Basic(U, Y, E, wnormtype, atol, vatol, rtol, vrtol, ignore_max, norm, norm_loc, norma, norma_loc, normr, normr_loc));
2502: PetscFunctionReturn(PETSC_SUCCESS);
2503: }