Actual source code: tagger.c
1: #include <petsc/private/vecimpl.h>
3: /*@
4: VecTaggerCreate - create a `VecTagger` context.
6: Collective
8: Input Parameter:
9: . comm - communicator on which the `VecTagger` will operate
11: Output Parameter:
12: . tagger - new Vec tagger context
14: Level: advanced
16: Notes:
17: This object is used to control the tagging/selection of index sets based on the values in a
18: vector. This is used, for example, in adaptive simulations when aspects are selected for
19: refinement or coarsening. The primary intent is that the selected index sets are based purely
20: on the values in the vector, though implementations that do not follow this intent are
21: possible.
23: Once a `VecTagger` is created (`VecTaggerCreate()`), optionally modified by options
24: (`VecTaggerSetFromOptions()`), and set up (`VecTaggerSetUp()`), it is applied to vectors with
25: `VecTaggerComputeIS()` to compute the selected index sets.
27: Provided implementations support tagging based on a box/interval of values
28: (`VECTAGGERABSOLUTE`), based on a box of values of relative to the range of values present in
29: the vector (`VECTAGGERRELATIVE`), based on where values fall in the cumulative distribution
30: of values in the vector (`VECTAGGERCDF`), and based on unions (`VECTAGGEROR`) or
31: intersections (`VECTAGGERAND`) of other criteria.
33: .seealso: `VecTagger`, `VecTaggerSetBlockSize()`, `VecTaggerSetFromOptions()`, `VecTaggerSetUp()`, `VecTaggerComputeIS()`, `VecTaggerComputeBoxes()`, `VecTaggerDestroy()`
34: @*/
35: PetscErrorCode VecTaggerCreate(MPI_Comm comm, VecTagger *tagger)
36: {
37: VecTagger b;
39: PetscFunctionBegin;
40: PetscAssertPointer(tagger, 2);
41: PetscCall(VecTaggerInitializePackage());
43: PetscCall(PetscHeaderCreate(b, VEC_TAGGER_CLASSID, "VecTagger", "Vec Tagger", "Vec", comm, VecTaggerDestroy, VecTaggerView));
44: b->blocksize = 1;
45: b->invert = PETSC_FALSE;
46: b->setupcalled = PETSC_FALSE;
47: *tagger = b;
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: /*@
52: VecTaggerSetType - set the Vec tagger implementation
54: Collective
56: Input Parameters:
57: + tagger - the `VecTagger` context
58: - type - a known method
60: Options Database Key:
61: . -vec_tagger_type <type> - Sets the method; use -help for a list
62: of available methods (for instance, absolute, relative, cdf, or, and)
64: Level: advanced
66: Notes:
67: See "include/petscvec.h" for available methods (for instance)
68: + `VECTAGGERABSOLUTE` - tag based on a box of values
69: . `VECTAGGERRELATIVE` - tag based on a box relative to the range of values present in the vector
70: . `VECTAGGERCDF` - tag based on a box in the cumulative distribution of values present in the vector
71: . `VECTAGGEROR` - tag based on the union of a set of `VecTagger` contexts
72: - `VECTAGGERAND` - tag based on the intersection of a set of other `VecTagger` contexts
74: .seealso: `VecTaggerType`, `VecTaggerCreate()`, `VecTagger`
75: @*/
76: PetscErrorCode VecTaggerSetType(VecTagger tagger, VecTaggerType type)
77: {
78: PetscBool match;
79: PetscErrorCode (*r)(VecTagger);
81: PetscFunctionBegin;
83: PetscAssertPointer(type, 2);
85: PetscCall(PetscObjectTypeCompare((PetscObject)tagger, type, &match));
86: if (match) PetscFunctionReturn(PETSC_SUCCESS);
88: PetscCall(PetscFunctionListFind(VecTaggerList, type, &r));
89: PetscCheck(r, PetscObjectComm((PetscObject)tagger), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unable to find requested VecTagger type %s", type);
90: /* Destroy the previous private VecTagger context */
91: PetscTryTypeMethod(tagger, destroy);
92: PetscCall(PetscMemzero(tagger->ops, sizeof(*tagger->ops)));
93: PetscCall(PetscObjectChangeTypeName((PetscObject)tagger, type));
94: tagger->ops->create = r;
95: PetscCall((*r)(tagger));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: /*@
100: VecTaggerGetType - Gets the `VecTaggerType` name (as a string) from the `VecTagger`.
102: Not Collective
104: Input Parameter:
105: . tagger - The `VecTagger` context
107: Output Parameter:
108: . type - The `VecTagger` type name
110: Level: advanced
112: .seealso: `VecTaggerSetType()`, `VecTaggerCreate()`, `VecTaggerSetFromOptions()`, `VecTagger`, `VecTaggerType`
113: @*/
114: PetscErrorCode VecTaggerGetType(VecTagger tagger, VecTaggerType *type)
115: {
116: PetscFunctionBegin;
118: PetscAssertPointer(type, 2);
119: PetscCall(VecTaggerRegisterAll());
120: *type = ((PetscObject)tagger)->type_name;
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
124: /*@
125: VecTaggerDestroy - destroy a `VecTagger` context
127: Collective
129: Input Parameter:
130: . tagger - address of tagger
132: Level: advanced
134: .seealso: `VecTaggerCreate()`, `VecTaggerSetType()`, `VecTagger`
135: @*/
136: PetscErrorCode VecTaggerDestroy(VecTagger *tagger)
137: {
138: PetscFunctionBegin;
139: if (!*tagger) PetscFunctionReturn(PETSC_SUCCESS);
141: if (--((PetscObject)*tagger)->refct > 0) {
142: *tagger = NULL;
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
145: PetscTryTypeMethod(*tagger, destroy);
146: PetscCall(PetscHeaderDestroy(tagger));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: /*@
151: VecTaggerSetUp - set up a `VecTagger` context
153: Collective
155: Input Parameter:
156: . tagger - Vec tagger object
158: Level: advanced
160: .seealso: `VecTaggerSetFromOptions()`, `VecTaggerSetType()`, `VecTagger`, `VecTaggerCreate()`
161: @*/
162: PetscErrorCode VecTaggerSetUp(VecTagger tagger)
163: {
164: PetscFunctionBegin;
165: if (tagger->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
166: if (!((PetscObject)tagger)->type_name) PetscCall(VecTaggerSetType(tagger, VECTAGGERABSOLUTE));
167: PetscTryTypeMethod(tagger, setup);
168: tagger->setupcalled = PETSC_TRUE;
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: /*@C
173: VecTaggerSetFromOptions - set `VecTagger` options using the options database
175: Logically Collective
177: Input Parameter:
178: . tagger - vec tagger
180: Options Database Keys:
181: + -vec_tagger_type - implementation type, see `VecTaggerSetType()`
182: . -vec_tagger_block_size - set the block size, see `VecTaggerSetBlockSize()`
183: - -vec_tagger_invert - invert the index set returned by `VecTaggerComputeIS()`
185: Level: advanced
187: .seealso: `VecTagger`, `VecTaggerCreate()`, `VecTaggerSetUp()`
189: @*/
190: PetscErrorCode VecTaggerSetFromOptions(VecTagger tagger)
191: {
192: VecTaggerType deft;
193: char type[256];
194: PetscBool flg;
196: PetscFunctionBegin;
198: PetscObjectOptionsBegin((PetscObject)tagger);
199: deft = ((PetscObject)tagger)->type_name ? ((PetscObject)tagger)->type_name : VECTAGGERABSOLUTE;
200: PetscCall(PetscOptionsFList("-vec_tagger_type", "VecTagger implementation type", "VecTaggerSetType", VecTaggerList, deft, type, 256, &flg));
201: PetscCall(VecTaggerSetType(tagger, flg ? type : deft));
202: PetscCall(PetscOptionsInt("-vec_tagger_block_size", "block size of the vectors the tagger operates on", "VecTaggerSetBlockSize", tagger->blocksize, &tagger->blocksize, NULL));
203: PetscCall(PetscOptionsBool("-vec_tagger_invert", "invert the set of indices returned by VecTaggerComputeIS()", "VecTaggerSetInvert", tagger->invert, &tagger->invert, NULL));
204: PetscTryTypeMethod(tagger, setfromoptions, PetscOptionsObject);
205: PetscOptionsEnd();
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: /*@
210: VecTaggerSetBlockSize - set the block size of the set of indices returned by `VecTaggerComputeIS()`.
212: Logically Collective
214: Input Parameters:
215: + tagger - vec tagger
216: - blocksize - block size of the criteria used to tagger vectors
218: Level: advanced
220: Notes:
221: Values greater than one are useful when there are multiple criteria for determining which
222: indices to include in the set. For example, consider adaptive mesh refinement in a
223: multiphysics problem, with metrics of solution quality for multiple fields measure on each
224: cell. The size of the vector will be `[numCells` * numFields]`; the `VecTagger` block size
225: should be `numFields`; `VecTaggerComputeIS()` will return indices in the range `[0,
226: numCells)`, i.e., one index is given for each block of values.
228: Note that the block size of the vector does not have to match this block size.
230: .seealso: `VecTaggerComputeIS()`, `VecTaggerGetBlockSize()`, `VecSetBlockSize()`, `VecGetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
231: @*/
232: PetscErrorCode VecTaggerSetBlockSize(VecTagger tagger, PetscInt blocksize)
233: {
234: PetscFunctionBegin;
237: tagger->blocksize = blocksize;
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: /*@
242: VecTaggerGetBlockSize - get the block size of the indices created by `VecTaggerComputeIS()`.
244: Logically Collective
246: Input Parameter:
247: . tagger - vec tagger
249: Output Parameter:
250: . blocksize - block size of the vectors the tagger operates on
252: Level: advanced
254: .seealso: `VecTaggerComputeIS()`, `VecTaggerSetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
255: @*/
256: PetscErrorCode VecTaggerGetBlockSize(VecTagger tagger, PetscInt *blocksize)
257: {
258: PetscFunctionBegin;
260: PetscAssertPointer(blocksize, 2);
261: *blocksize = tagger->blocksize;
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: /*@
266: VecTaggerSetInvert - If the tagged index sets are based on boxes that can be returned by `VecTaggerComputeBoxes()`,
267: then this option inverts values used to compute the IS, i.e., from being in the union of the boxes to being in the
268: intersection of their exteriors.
270: Logically Collective
272: Input Parameters:
273: + tagger - vec tagger
274: - invert - `PETSC_TRUE` to invert, `PETSC_FALSE` to use the indices as is
276: Level: advanced
278: .seealso: `VecTaggerComputeIS()`, `VecTaggerGetInvert()`, `VecTagger`, `VecTaggerCreate()`
279: @*/
280: PetscErrorCode VecTaggerSetInvert(VecTagger tagger, PetscBool invert)
281: {
282: PetscFunctionBegin;
285: tagger->invert = invert;
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: /*@
290: VecTaggerGetInvert - get whether the set of indices returned by `VecTaggerComputeIS()` are inverted
292: Logically Collective
294: Input Parameter:
295: . tagger - vec tagger
297: Output Parameter:
298: . invert - `PETSC_TRUE` to invert, `PETSC_FALSE` to use the indices as is
300: Level: advanced
302: .seealso: `VecTaggerComputeIS()`, `VecTaggerSetInvert()`, `VecTagger`, `VecTaggerCreate()`
303: @*/
304: PetscErrorCode VecTaggerGetInvert(VecTagger tagger, PetscBool *invert)
305: {
306: PetscFunctionBegin;
308: PetscAssertPointer(invert, 2);
309: *invert = tagger->invert;
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: /*@
314: VecTaggerView - view a `VecTagger` context
316: Collective
318: Input Parameters:
319: + tagger - vec tagger
320: - viewer - viewer to display tagger, for example `PETSC_VIEWER_STDOUT_WORLD`
322: Level: advanced
324: .seealso: `VecTaggerCreate()`, `VecTagger`
325: @*/
326: PetscErrorCode VecTaggerView(VecTagger tagger, PetscViewer viewer)
327: {
328: PetscBool iascii;
330: PetscFunctionBegin;
332: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tagger), &viewer));
334: PetscCheckSameComm(tagger, 1, viewer, 2);
335: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
336: if (iascii) {
337: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tagger, viewer));
338: PetscCall(PetscViewerASCIIPushTab(viewer));
339: PetscCall(PetscViewerASCIIPrintf(viewer, "Block size: %" PetscInt_FMT "\n", tagger->blocksize));
340: PetscTryTypeMethod(tagger, view, viewer);
341: if (tagger->invert) PetscCall(PetscViewerASCIIPrintf(viewer, "Inverting ISs.\n"));
342: PetscCall(PetscViewerASCIIPopTab(viewer));
343: }
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: /*@C
348: VecTaggerComputeBoxes - If the tagged index set can be summarized as a list of boxes of values, returns that list, otherwise returns
349: in listed `PETSC_FALSE`
351: Collective
353: Input Parameters:
354: + tagger - the `VecTagger` context
355: - vec - the vec to tag
357: Output Parameters:
358: + numBoxes - the number of boxes in the tag definition
359: . boxes - a newly allocated list of boxes. This is a flat array of (BlockSize * `numBoxe`s) pairs that the user can free with `PetscFree()`.
360: - listed - `PETSC_TRUE` if a list was created, pass in `NULL` if not needed
362: Level: advanced
364: Note:
365: A value is tagged if it is in any of the boxes, unless the tagger has been inverted (see `VecTaggerSetInvert()`/`VecTaggerGetInvert()`), in which case a value is tagged if it is in none of the boxes.
367: .seealso: `VecTaggerComputeIS()`, `VecTagger`, `VecTaggerCreate()`
368: @*/
369: PetscErrorCode VecTaggerComputeBoxes(VecTagger tagger, Vec vec, PetscInt *numBoxes, VecTaggerBox *boxes[], PetscBool *listed)
370: {
371: PetscInt vls, tbs;
373: PetscFunctionBegin;
376: PetscAssertPointer(numBoxes, 3);
377: PetscAssertPointer(boxes, 4);
378: PetscCall(VecGetLocalSize(vec, &vls));
379: PetscCall(VecTaggerGetBlockSize(tagger, &tbs));
380: PetscCheck(vls % tbs == 0, PetscObjectComm((PetscObject)tagger), PETSC_ERR_ARG_INCOMP, "vec local size %" PetscInt_FMT " is not a multiple of tagger block size %" PetscInt_FMT, vls, tbs);
381: if (tagger->ops->computeboxes) {
382: *listed = PETSC_TRUE;
383: PetscUseTypeMethod(tagger, computeboxes, vec, numBoxes, boxes, listed);
384: } else *listed = PETSC_FALSE;
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: /*@C
389: VecTaggerComputeIS - Use a `VecTagger` context to tag a set of indices based on a vector's values
391: Collective
393: Input Parameters:
394: + tagger - the `VecTagger` context
395: - vec - the vec to tag
397: Output Parameters:
398: + is - a list of the indices tagged by the tagger, i.e., if the number of local indices will be n / bs, where n is `VecGetLocalSize()` and bs is `VecTaggerGetBlockSize()`.
399: - listed - routine was able to compute the `IS`, pass in `NULL` if not needed
401: Level: advanced
403: .seealso: `VecTaggerComputeBoxes()`, `VecTagger`, `VecTaggerCreate()`
404: @*/
405: PetscErrorCode VecTaggerComputeIS(VecTagger tagger, Vec vec, IS is[], PetscBool *listed)
406: {
407: PetscInt vls, tbs;
409: PetscFunctionBegin;
412: PetscAssertPointer(is, 3);
413: PetscCall(VecGetLocalSize(vec, &vls));
414: PetscCall(VecTaggerGetBlockSize(tagger, &tbs));
415: PetscCheck(vls % tbs == 0, PetscObjectComm((PetscObject)tagger), PETSC_ERR_ARG_INCOMP, "vec local size %" PetscInt_FMT " is not a multiple of tagger block size %" PetscInt_FMT, vls, tbs);
416: if (tagger->ops->computeis) {
417: PetscUseTypeMethod(tagger, computeis, vec, is, listed);
418: } else if (listed) *listed = PETSC_FALSE;
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: PetscErrorCode VecTaggerComputeIS_FromBoxes(VecTagger tagger, Vec vec, IS *is, PetscBool *listed)
423: {
424: PetscInt numBoxes;
425: VecTaggerBox *boxes;
426: PetscInt numTagged, offset;
427: PetscInt *tagged;
428: PetscInt bs, b, i, j, k, n;
429: PetscBool invert;
430: const PetscScalar *vecArray;
431: PetscBool boxlisted;
433: PetscFunctionBegin;
434: PetscCall(VecTaggerGetBlockSize(tagger, &bs));
435: PetscCall(VecTaggerComputeBoxes(tagger, vec, &numBoxes, &boxes, &boxlisted));
436: if (!boxlisted) {
437: if (listed) *listed = PETSC_FALSE;
438: PetscFunctionReturn(PETSC_SUCCESS);
439: }
440: PetscCall(VecGetArrayRead(vec, &vecArray));
441: PetscCall(VecGetLocalSize(vec, &n));
442: invert = tagger->invert;
443: numTagged = 0;
444: offset = 0;
445: tagged = NULL;
446: PetscCheck(n % bs == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "blocksize %" PetscInt_FMT " does not divide vector length %" PetscInt_FMT, bs, n);
447: n /= bs;
448: for (i = 0; i < 2; i++) {
449: if (i) PetscCall(PetscMalloc1(numTagged, &tagged));
450: for (j = 0; j < n; j++) {
451: for (k = 0; k < numBoxes; k++) {
452: for (b = 0; b < bs; b++) {
453: PetscScalar val = vecArray[j * bs + b];
454: PetscInt l = k * bs + b;
455: VecTaggerBox box;
456: PetscBool in;
458: box = boxes[l];
459: #if !defined(PETSC_USE_COMPLEX)
460: in = (PetscBool)((box.min <= val) && (val <= box.max));
461: #else
462: in = (PetscBool)((PetscRealPart(box.min) <= PetscRealPart(val)) && (PetscImaginaryPart(box.min) <= PetscImaginaryPart(val)) && (PetscRealPart(val) <= PetscRealPart(box.max)) && (PetscImaginaryPart(val) <= PetscImaginaryPart(box.max)));
463: #endif
464: if (!in) break;
465: }
466: if (b == bs) break;
467: }
468: if ((PetscBool)(k < numBoxes) ^ invert) {
469: if (!i) numTagged++;
470: else tagged[offset++] = j;
471: }
472: }
473: }
474: PetscCall(VecRestoreArrayRead(vec, &vecArray));
475: PetscCall(PetscFree(boxes));
476: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)vec), numTagged, tagged, PETSC_OWN_POINTER, is));
477: PetscCall(ISSort(*is));
478: if (listed) *listed = PETSC_TRUE;
479: PetscFunctionReturn(PETSC_SUCCESS);
480: }