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));
45: b->blocksize = 1;
46: b->invert = PETSC_FALSE;
47: b->setupcalled = PETSC_FALSE;
49: *tagger = b;
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: /*@
54: VecTaggerSetType - set the Vec tagger implementation
56: Collective
58: Input Parameters:
59: + tagger - the `VecTagger` context
60: - type - a known method
62: Options Database Key:
63: . -vec_tagger_type <type> - Sets the method; use -help for a list
64: of available methods (for instance, absolute, relative, cdf, or, and)
66: Level: advanced
68: Notes:
69: See "include/petscvec.h" for available methods (for instance)
70: + `VECTAGGERABSOLUTE` - tag based on a box of values
71: . `VECTAGGERRELATIVE` - tag based on a box relative to the range of values present in the vector
72: . `VECTAGGERCDF` - tag based on a box in the cumulative distribution of values present in the vector
73: . `VECTAGGEROR` - tag based on the union of a set of `VecTagger` contexts
74: - `VECTAGGERAND` - tag based on the intersection of a set of other `VecTagger` contexts
76: .seealso: `VecTaggerType`, `VecTaggerCreate()`, `VecTagger`
77: @*/
78: PetscErrorCode VecTaggerSetType(VecTagger tagger, VecTaggerType type)
79: {
80: PetscBool match;
81: PetscErrorCode (*r)(VecTagger);
83: PetscFunctionBegin;
85: PetscAssertPointer(type, 2);
87: PetscCall(PetscObjectTypeCompare((PetscObject)tagger, type, &match));
88: if (match) PetscFunctionReturn(PETSC_SUCCESS);
90: PetscCall(PetscFunctionListFind(VecTaggerList, type, &r));
91: PetscCheck(r, PetscObjectComm((PetscObject)tagger), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unable to find requested VecTagger type %s", type);
92: /* Destroy the previous private VecTagger context */
93: PetscTryTypeMethod(tagger, destroy);
94: PetscCall(PetscMemzero(tagger->ops, sizeof(*tagger->ops)));
95: PetscCall(PetscObjectChangeTypeName((PetscObject)tagger, type));
96: tagger->ops->create = r;
97: PetscCall((*r)(tagger));
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: /*@
102: VecTaggerGetType - Gets the `VecTaggerType` name (as a string) from the `VecTagger`.
104: Not Collective
106: Input Parameter:
107: . tagger - The `VecTagger` context
109: Output Parameter:
110: . type - The `VecTagger` type name
112: Level: advanced
114: .seealso: `VecTaggerSetType()`, `VecTaggerCreate()`, `VecTaggerSetFromOptions()`, `VecTagger`, `VecTaggerType`
115: @*/
116: PetscErrorCode VecTaggerGetType(VecTagger tagger, VecTaggerType *type)
117: {
118: PetscFunctionBegin;
120: PetscAssertPointer(type, 2);
121: PetscCall(VecTaggerRegisterAll());
122: *type = ((PetscObject)tagger)->type_name;
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: /*@
127: VecTaggerDestroy - destroy a `VecTagger` context
129: Collective
131: Input Parameter:
132: . tagger - address of tagger
134: Level: advanced
136: .seealso: `VecTaggerCreate()`, `VecTaggerSetType()`, `VecTagger`
137: @*/
138: PetscErrorCode VecTaggerDestroy(VecTagger *tagger)
139: {
140: PetscFunctionBegin;
141: if (!*tagger) PetscFunctionReturn(PETSC_SUCCESS);
143: if (--((PetscObject)*tagger)->refct > 0) {
144: *tagger = NULL;
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
147: PetscTryTypeMethod(*tagger, destroy);
148: PetscCall(PetscHeaderDestroy(tagger));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: /*@
153: VecTaggerSetUp - set up a `VecTagger` context
155: Collective
157: Input Parameter:
158: . tagger - Vec tagger object
160: Level: advanced
162: .seealso: `VecTaggerSetFromOptions()`, `VecTaggerSetType()`, `VecTagger`, `VecTaggerCreate()`
163: @*/
164: PetscErrorCode VecTaggerSetUp(VecTagger tagger)
165: {
166: PetscFunctionBegin;
167: if (tagger->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
168: if (!((PetscObject)tagger)->type_name) PetscCall(VecTaggerSetType(tagger, VECTAGGERABSOLUTE));
169: PetscTryTypeMethod(tagger, setup);
170: tagger->setupcalled = PETSC_TRUE;
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: /*@C
175: VecTaggerSetFromOptions - set `VecTagger` options using the options database
177: Logically Collective
179: Input Parameter:
180: . tagger - vec tagger
182: Options Database Keys:
183: + -vec_tagger_type - implementation type, see `VecTaggerSetType()`
184: . -vec_tagger_block_size - set the block size, see `VecTaggerSetBlockSize()`
185: - -vec_tagger_invert - invert the index set returned by `VecTaggerComputeIS()`
187: Level: advanced
189: .seealso: `VecTagger`, `VecTaggerCreate()`, `VecTaggerSetUp()`
191: @*/
192: PetscErrorCode VecTaggerSetFromOptions(VecTagger tagger)
193: {
194: VecTaggerType deft;
195: char type[256];
196: PetscBool flg;
198: PetscFunctionBegin;
200: PetscObjectOptionsBegin((PetscObject)tagger);
201: deft = ((PetscObject)tagger)->type_name ? ((PetscObject)tagger)->type_name : VECTAGGERABSOLUTE;
202: PetscCall(PetscOptionsFList("-vec_tagger_type", "VecTagger implementation type", "VecTaggerSetType", VecTaggerList, deft, type, 256, &flg));
203: PetscCall(VecTaggerSetType(tagger, flg ? type : deft));
204: PetscCall(PetscOptionsInt("-vec_tagger_block_size", "block size of the vectors the tagger operates on", "VecTaggerSetBlockSize", tagger->blocksize, &tagger->blocksize, NULL));
205: PetscCall(PetscOptionsBool("-vec_tagger_invert", "invert the set of indices returned by VecTaggerComputeIS()", "VecTaggerSetInvert", tagger->invert, &tagger->invert, NULL));
206: PetscTryTypeMethod(tagger, setfromoptions, PetscOptionsObject);
207: PetscOptionsEnd();
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: /*@
212: VecTaggerSetBlockSize - set the block size of the set of indices returned by `VecTaggerComputeIS()`.
214: Logically Collective
216: Input Parameters:
217: + tagger - vec tagger
218: - blocksize - block size of the criteria used to tagger vectors
220: Level: advanced
222: Notes:
223: Values greater than one are useful when there are multiple criteria for determining which
224: indices to include in the set. For example, consider adaptive mesh refinement in a
225: multiphysics problem, with metrics of solution quality for multiple fields measure on each
226: cell. The size of the vector will be `[numCells` * numFields]`; the `VecTagger` block size
227: should be `numFields`; `VecTaggerComputeIS()` will return indices in the range `[0,
228: numCells)`, i.e., one index is given for each block of values.
230: Note that the block size of the vector does not have to match this block size.
232: .seealso: `VecTaggerComputeIS()`, `VecTaggerGetBlockSize()`, `VecSetBlockSize()`, `VecGetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
233: @*/
234: PetscErrorCode VecTaggerSetBlockSize(VecTagger tagger, PetscInt blocksize)
235: {
236: PetscFunctionBegin;
239: tagger->blocksize = blocksize;
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: /*@
244: VecTaggerGetBlockSize - get the block size of the indices created by `VecTaggerComputeIS()`.
246: Logically Collective
248: Input Parameter:
249: . tagger - vec tagger
251: Output Parameter:
252: . blocksize - block size of the vectors the tagger operates on
254: Level: advanced
256: .seealso: `VecTaggerComputeIS()`, `VecTaggerSetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
257: @*/
258: PetscErrorCode VecTaggerGetBlockSize(VecTagger tagger, PetscInt *blocksize)
259: {
260: PetscFunctionBegin;
262: PetscAssertPointer(blocksize, 2);
263: *blocksize = tagger->blocksize;
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: /*@
268: VecTaggerSetInvert - If the tagged index sets are based on boxes that can be returned by `VecTaggerComputeBoxes()`,
269: then this option inverts values used to compute the IS, i.e., from being in the union of the boxes to being in the
270: intersection of their exteriors.
272: Logically Collective
274: Input Parameters:
275: + tagger - vec tagger
276: - invert - `PETSC_TRUE` to invert, `PETSC_FALSE` to use the indices as is
278: Level: advanced
280: .seealso: `VecTaggerComputeIS()`, `VecTaggerGetInvert()`, `VecTagger`, `VecTaggerCreate()`
281: @*/
282: PetscErrorCode VecTaggerSetInvert(VecTagger tagger, PetscBool invert)
283: {
284: PetscFunctionBegin;
287: tagger->invert = invert;
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: /*@
292: VecTaggerGetInvert - get whether the set of indices returned by `VecTaggerComputeIS()` are inverted
294: Logically Collective
296: Input Parameter:
297: . tagger - vec tagger
299: Output Parameter:
300: . invert - `PETSC_TRUE` to invert, `PETSC_FALSE` to use the indices as is
302: Level: advanced
304: .seealso: `VecTaggerComputeIS()`, `VecTaggerSetInvert()`, `VecTagger`, `VecTaggerCreate()`
305: @*/
306: PetscErrorCode VecTaggerGetInvert(VecTagger tagger, PetscBool *invert)
307: {
308: PetscFunctionBegin;
310: PetscAssertPointer(invert, 2);
311: *invert = tagger->invert;
312: PetscFunctionReturn(PETSC_SUCCESS);
313: }
315: /*@C
316: VecTaggerView - view a `VecTagger` context
318: Collective
320: Input Parameters:
321: + tagger - vec tagger
322: - viewer - viewer to display tagger, for example `PETSC_VIEWER_STDOUT_WORLD`
324: Level: advanced
326: .seealso: `VecTaggerCreate()`, `VecTagger`
327: @*/
328: PetscErrorCode VecTaggerView(VecTagger tagger, PetscViewer viewer)
329: {
330: PetscBool iascii;
332: PetscFunctionBegin;
334: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tagger), &viewer));
336: PetscCheckSameComm(tagger, 1, viewer, 2);
337: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
338: if (iascii) {
339: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tagger, viewer));
340: PetscCall(PetscViewerASCIIPushTab(viewer));
341: PetscCall(PetscViewerASCIIPrintf(viewer, "Block size: %" PetscInt_FMT "\n", tagger->blocksize));
342: PetscTryTypeMethod(tagger, view, viewer);
343: if (tagger->invert) PetscCall(PetscViewerASCIIPrintf(viewer, "Inverting ISs.\n"));
344: PetscCall(PetscViewerASCIIPopTab(viewer));
345: }
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: /*@C
350: VecTaggerComputeBoxes - If the tagged index set can be summarized as a list of boxes of values, returns that list, otherwise returns
351: in listed `PETSC_FALSE`
353: Collective
355: Input Parameters:
356: + tagger - the `VecTagger` context
357: - vec - the vec to tag
359: Output Parameters:
360: + numBoxes - the number of boxes in the tag definition
361: . boxes - a newly allocated list of boxes. This is a flat array of (BlockSize * `numBoxe`s) pairs that the user can free with `PetscFree()`.
362: - listed - `PETSC_TRUE` if a list was created, pass in `NULL` if not needed
364: Level: advanced
366: Note:
367: 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.
369: .seealso: `VecTaggerComputeIS()`, `VecTagger`, `VecTaggerCreate()`
370: @*/
371: PetscErrorCode VecTaggerComputeBoxes(VecTagger tagger, Vec vec, PetscInt *numBoxes, VecTaggerBox *boxes[], PetscBool *listed)
372: {
373: PetscInt vls, tbs;
375: PetscFunctionBegin;
378: PetscAssertPointer(numBoxes, 3);
379: PetscAssertPointer(boxes, 4);
380: PetscCall(VecGetLocalSize(vec, &vls));
381: PetscCall(VecTaggerGetBlockSize(tagger, &tbs));
382: 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);
383: if (tagger->ops->computeboxes) {
384: *listed = PETSC_TRUE;
385: PetscUseTypeMethod(tagger, computeboxes, vec, numBoxes, boxes, listed);
386: } else *listed = PETSC_FALSE;
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: /*@C
391: VecTaggerComputeIS - Use a `VecTagger` context to tag a set of indices based on a vector's values
393: Collective
395: Input Parameters:
396: + tagger - the `VecTagger` context
397: - vec - the vec to tag
399: Output Parameters:
400: + 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()`.
401: - listed - routine was able to compute the `IS`, pass in `NULL` if not needed
403: Level: advanced
405: .seealso: `VecTaggerComputeBoxes()`, `VecTagger`, `VecTaggerCreate()`
406: @*/
407: PetscErrorCode VecTaggerComputeIS(VecTagger tagger, Vec vec, IS is[], PetscBool *listed)
408: {
409: PetscInt vls, tbs;
411: PetscFunctionBegin;
414: PetscAssertPointer(is, 3);
415: PetscCall(VecGetLocalSize(vec, &vls));
416: PetscCall(VecTaggerGetBlockSize(tagger, &tbs));
417: 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);
418: if (tagger->ops->computeis) {
419: PetscUseTypeMethod(tagger, computeis, vec, is, listed);
420: } else if (listed) *listed = PETSC_FALSE;
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: PetscErrorCode VecTaggerComputeIS_FromBoxes(VecTagger tagger, Vec vec, IS *is, PetscBool *listed)
425: {
426: PetscInt numBoxes;
427: VecTaggerBox *boxes;
428: PetscInt numTagged, offset;
429: PetscInt *tagged;
430: PetscInt bs, b, i, j, k, n;
431: PetscBool invert;
432: const PetscScalar *vecArray;
433: PetscBool boxlisted;
435: PetscFunctionBegin;
436: PetscCall(VecTaggerGetBlockSize(tagger, &bs));
437: PetscCall(VecTaggerComputeBoxes(tagger, vec, &numBoxes, &boxes, &boxlisted));
438: if (!boxlisted) {
439: if (listed) *listed = PETSC_FALSE;
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
442: PetscCall(VecGetArrayRead(vec, &vecArray));
443: PetscCall(VecGetLocalSize(vec, &n));
444: invert = tagger->invert;
445: numTagged = 0;
446: offset = 0;
447: tagged = NULL;
448: PetscCheck(n % bs == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "blocksize %" PetscInt_FMT " does not divide vector length %" PetscInt_FMT, bs, n);
449: n /= bs;
450: for (i = 0; i < 2; i++) {
451: if (i) PetscCall(PetscMalloc1(numTagged, &tagged));
452: for (j = 0; j < n; j++) {
453: for (k = 0; k < numBoxes; k++) {
454: for (b = 0; b < bs; b++) {
455: PetscScalar val = vecArray[j * bs + b];
456: PetscInt l = k * bs + b;
457: VecTaggerBox box;
458: PetscBool in;
460: box = boxes[l];
461: #if !defined(PETSC_USE_COMPLEX)
462: in = (PetscBool)((box.min <= val) && (val <= box.max));
463: #else
464: in = (PetscBool)((PetscRealPart(box.min) <= PetscRealPart(val)) && (PetscImaginaryPart(box.min) <= PetscImaginaryPart(val)) && (PetscRealPart(val) <= PetscRealPart(box.max)) && (PetscImaginaryPart(val) <= PetscImaginaryPart(box.max)));
465: #endif
466: if (!in) break;
467: }
468: if (b == bs) break;
469: }
470: if ((PetscBool)(k < numBoxes) ^ invert) {
471: if (!i) numTagged++;
472: else tagged[offset++] = j;
473: }
474: }
475: }
476: PetscCall(VecRestoreArrayRead(vec, &vecArray));
477: PetscCall(PetscFree(boxes));
478: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)vec), numTagged, tagged, PETSC_OWN_POINTER, is));
479: PetscCall(ISSort(*is));
480: if (listed) *listed = PETSC_TRUE;
481: PetscFunctionReturn(PETSC_SUCCESS);
482: }