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; see `VecTaggerType`
63: Level: advanced
65: Notes:
66: See "include/petscvec.h" for available methods (for instance)
67: + `VECTAGGERABSOLUTE` - tag based on a box of values
68: . `VECTAGGERRELATIVE` - tag based on a box relative to the range of values present in the vector
69: . `VECTAGGERCDF` - tag based on a box in the cumulative distribution of values present in the vector
70: . `VECTAGGEROR` - tag based on the union of a set of `VecTagger` contexts
71: - `VECTAGGERAND` - tag based on the intersection of a set of other `VecTagger` contexts
73: .seealso: `VecTaggerType`, `VecTaggerCreate()`, `VecTagger`
74: @*/
75: PetscErrorCode VecTaggerSetType(VecTagger tagger, VecTaggerType type)
76: {
77: PetscBool match;
78: PetscErrorCode (*r)(VecTagger);
80: PetscFunctionBegin;
82: PetscAssertPointer(type, 2);
84: PetscCall(PetscObjectTypeCompare((PetscObject)tagger, type, &match));
85: if (match) PetscFunctionReturn(PETSC_SUCCESS);
87: PetscCall(PetscFunctionListFind(VecTaggerList, type, &r));
88: PetscCheck(r, PetscObjectComm((PetscObject)tagger), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unable to find requested VecTagger type %s", type);
89: /* Destroy the previous private VecTagger context */
90: PetscTryTypeMethod(tagger, destroy);
91: PetscCall(PetscMemzero(tagger->ops, sizeof(*tagger->ops)));
92: PetscCall(PetscObjectChangeTypeName((PetscObject)tagger, type));
93: tagger->ops->create = r;
94: PetscCall((*r)(tagger));
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: /*@
99: VecTaggerGetType - Gets the `VecTaggerType` name (as a string) from the `VecTagger`.
101: Not Collective
103: Input Parameter:
104: . tagger - The `VecTagger` context
106: Output Parameter:
107: . type - The `VecTagger` type name
109: Level: advanced
111: .seealso: `VecTaggerSetType()`, `VecTaggerCreate()`, `VecTaggerSetFromOptions()`, `VecTagger`, `VecTaggerType`
112: @*/
113: PetscErrorCode VecTaggerGetType(VecTagger tagger, VecTaggerType *type)
114: {
115: PetscFunctionBegin;
117: PetscAssertPointer(type, 2);
118: PetscCall(VecTaggerRegisterAll());
119: *type = ((PetscObject)tagger)->type_name;
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: /*@
124: VecTaggerDestroy - destroy a `VecTagger` context
126: Collective
128: Input Parameter:
129: . tagger - address of tagger
131: Level: advanced
133: .seealso: `VecTaggerCreate()`, `VecTaggerSetType()`, `VecTagger`
134: @*/
135: PetscErrorCode VecTaggerDestroy(VecTagger *tagger)
136: {
137: PetscFunctionBegin;
138: if (!*tagger) PetscFunctionReturn(PETSC_SUCCESS);
140: if (--((PetscObject)*tagger)->refct > 0) {
141: *tagger = NULL;
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
144: PetscTryTypeMethod(*tagger, destroy);
145: PetscCall(PetscHeaderDestroy(tagger));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*@
150: VecTaggerSetUp - set up a `VecTagger` context
152: Collective
154: Input Parameter:
155: . tagger - Vec tagger object
157: Level: advanced
159: .seealso: `VecTaggerSetFromOptions()`, `VecTaggerSetType()`, `VecTagger`, `VecTaggerCreate()`
160: @*/
161: PetscErrorCode VecTaggerSetUp(VecTagger tagger)
162: {
163: PetscFunctionBegin;
164: if (tagger->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
165: if (!((PetscObject)tagger)->type_name) PetscCall(VecTaggerSetType(tagger, VECTAGGERABSOLUTE));
166: PetscTryTypeMethod(tagger, setup);
167: tagger->setupcalled = PETSC_TRUE;
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: /*@C
172: VecTaggerSetFromOptions - set `VecTagger` options using the options database
174: Logically Collective
176: Input Parameter:
177: . tagger - vec tagger
179: Options Database Keys:
180: + -vec_tagger_type - implementation type, see `VecTaggerSetType()`
181: . -vec_tagger_block_size - set the block size, see `VecTaggerSetBlockSize()`
182: - -vec_tagger_invert - invert the index set returned by `VecTaggerComputeIS()`
184: Level: advanced
186: .seealso: `VecTagger`, `VecTaggerCreate()`, `VecTaggerSetUp()`
187: @*/
188: PetscErrorCode VecTaggerSetFromOptions(VecTagger tagger)
189: {
190: VecTaggerType deft;
191: char type[256];
192: PetscBool flg;
194: PetscFunctionBegin;
196: PetscObjectOptionsBegin((PetscObject)tagger);
197: deft = ((PetscObject)tagger)->type_name ? ((PetscObject)tagger)->type_name : VECTAGGERABSOLUTE;
198: PetscCall(PetscOptionsFList("-vec_tagger_type", "VecTagger implementation type", "VecTaggerSetType", VecTaggerList, deft, type, 256, &flg));
199: PetscCall(VecTaggerSetType(tagger, flg ? type : deft));
200: PetscCall(PetscOptionsInt("-vec_tagger_block_size", "block size of the vectors the tagger operates on", "VecTaggerSetBlockSize", tagger->blocksize, &tagger->blocksize, NULL));
201: PetscCall(PetscOptionsBool("-vec_tagger_invert", "invert the set of indices returned by VecTaggerComputeIS()", "VecTaggerSetInvert", tagger->invert, &tagger->invert, NULL));
202: PetscTryTypeMethod(tagger, setfromoptions, PetscOptionsObject);
203: PetscOptionsEnd();
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: /*@
208: VecTaggerSetBlockSize - set the block size of the set of indices returned by `VecTaggerComputeIS()`.
210: Logically Collective
212: Input Parameters:
213: + tagger - vec tagger
214: - blocksize - block size of the criteria used to tagger vectors
216: Level: advanced
218: Notes:
219: Values greater than one are useful when there are multiple criteria for determining which
220: indices to include in the set. For example, consider adaptive mesh refinement in a
221: multiphysics problem, with metrics of solution quality for multiple fields measure on each
222: cell. The size of the vector will be `[numCells` * numFields]`; the `VecTagger` block size
223: should be `numFields`; `VecTaggerComputeIS()` will return indices in the range `[0,
224: numCells)`, i.e., one index is given for each block of values.
226: Note that the block size of the vector does not have to match this block size.
228: .seealso: `VecTaggerComputeIS()`, `VecTaggerGetBlockSize()`, `VecSetBlockSize()`, `VecGetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
229: @*/
230: PetscErrorCode VecTaggerSetBlockSize(VecTagger tagger, PetscInt blocksize)
231: {
232: PetscFunctionBegin;
235: tagger->blocksize = blocksize;
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: /*@
240: VecTaggerGetBlockSize - get the block size of the indices created by `VecTaggerComputeIS()`.
242: Logically Collective
244: Input Parameter:
245: . tagger - vec tagger
247: Output Parameter:
248: . blocksize - block size of the vectors the tagger operates on
250: Level: advanced
252: .seealso: `VecTaggerComputeIS()`, `VecTaggerSetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
253: @*/
254: PetscErrorCode VecTaggerGetBlockSize(VecTagger tagger, PetscInt *blocksize)
255: {
256: PetscFunctionBegin;
258: PetscAssertPointer(blocksize, 2);
259: *blocksize = tagger->blocksize;
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: /*@
264: VecTaggerSetInvert - If the tagged index sets are based on boxes that can be returned by `VecTaggerComputeBoxes()`,
265: then this option inverts values used to compute the IS, i.e., from being in the union of the boxes to being in the
266: intersection of their exteriors.
268: Logically Collective
270: Input Parameters:
271: + tagger - vec tagger
272: - invert - `PETSC_TRUE` to invert, `PETSC_FALSE` to use the indices as is
274: Level: advanced
276: .seealso: `VecTaggerComputeIS()`, `VecTaggerGetInvert()`, `VecTagger`, `VecTaggerCreate()`
277: @*/
278: PetscErrorCode VecTaggerSetInvert(VecTagger tagger, PetscBool invert)
279: {
280: PetscFunctionBegin;
283: tagger->invert = invert;
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: /*@
288: VecTaggerGetInvert - get whether the set of indices returned by `VecTaggerComputeIS()` are inverted
290: Logically Collective
292: Input Parameter:
293: . tagger - vec tagger
295: Output Parameter:
296: . invert - `PETSC_TRUE` to invert, `PETSC_FALSE` to use the indices as is
298: Level: advanced
300: .seealso: `VecTaggerComputeIS()`, `VecTaggerSetInvert()`, `VecTagger`, `VecTaggerCreate()`
301: @*/
302: PetscErrorCode VecTaggerGetInvert(VecTagger tagger, PetscBool *invert)
303: {
304: PetscFunctionBegin;
306: PetscAssertPointer(invert, 2);
307: *invert = tagger->invert;
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: /*@
312: VecTaggerView - view a `VecTagger` context
314: Collective
316: Input Parameters:
317: + tagger - vec tagger
318: - viewer - viewer to display tagger, for example `PETSC_VIEWER_STDOUT_WORLD`
320: Level: advanced
322: .seealso: `VecTaggerCreate()`, `VecTagger`
323: @*/
324: PetscErrorCode VecTaggerView(VecTagger tagger, PetscViewer viewer)
325: {
326: PetscBool isascii;
328: PetscFunctionBegin;
330: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tagger), &viewer));
332: PetscCheckSameComm(tagger, 1, viewer, 2);
333: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
334: if (isascii) {
335: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tagger, viewer));
336: PetscCall(PetscViewerASCIIPushTab(viewer));
337: PetscCall(PetscViewerASCIIPrintf(viewer, "Block size: %" PetscInt_FMT "\n", tagger->blocksize));
338: PetscTryTypeMethod(tagger, view, viewer);
339: if (tagger->invert) PetscCall(PetscViewerASCIIPrintf(viewer, "Inverting ISs.\n"));
340: PetscCall(PetscViewerASCIIPopTab(viewer));
341: }
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: /*@C
346: VecTaggerComputeBoxes - If the tagged index set can be summarized as a list of boxes of values, returns that list, otherwise returns
347: in listed `PETSC_FALSE`
349: Collective
351: Input Parameters:
352: + tagger - the `VecTagger` context
353: - vec - the vec to tag
355: Output Parameters:
356: + numBoxes - the number of boxes in the tag definition
357: . boxes - a newly allocated list of boxes. This is a flat array of (BlockSize * `numBoxe`s) pairs that the user can free with `PetscFree()`.
358: - listed - `PETSC_TRUE` if a list was created, pass in `NULL` if not needed
360: Level: advanced
362: Note:
363: 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.
365: .seealso: `VecTaggerComputeIS()`, `VecTagger`, `VecTaggerCreate()`
366: @*/
367: PetscErrorCode VecTaggerComputeBoxes(VecTagger tagger, Vec vec, PetscInt *numBoxes, VecTaggerBox *boxes[], PetscBool *listed)
368: {
369: PetscInt vls, tbs;
371: PetscFunctionBegin;
374: PetscAssertPointer(numBoxes, 3);
375: PetscAssertPointer(boxes, 4);
376: PetscCall(VecGetLocalSize(vec, &vls));
377: PetscCall(VecTaggerGetBlockSize(tagger, &tbs));
378: 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);
379: if (tagger->ops->computeboxes) {
380: *listed = PETSC_TRUE;
381: PetscUseTypeMethod(tagger, computeboxes, vec, numBoxes, boxes, listed);
382: } else *listed = PETSC_FALSE;
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: /*@C
387: VecTaggerComputeIS - Use a `VecTagger` context to tag a set of indices based on a vector's values
389: Collective
391: Input Parameters:
392: + tagger - the `VecTagger` context
393: - vec - the vec to tag
395: Output Parameters:
396: + 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()`.
397: - listed - routine was able to compute the `IS`, pass in `NULL` if not needed
399: Level: advanced
401: .seealso: `VecTaggerComputeBoxes()`, `VecTagger`, `VecTaggerCreate()`
402: @*/
403: PetscErrorCode VecTaggerComputeIS(VecTagger tagger, Vec vec, IS is[], PetscBool *listed)
404: {
405: PetscInt vls, tbs;
407: PetscFunctionBegin;
410: PetscAssertPointer(is, 3);
411: PetscCall(VecGetLocalSize(vec, &vls));
412: PetscCall(VecTaggerGetBlockSize(tagger, &tbs));
413: 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);
414: if (tagger->ops->computeis) {
415: PetscUseTypeMethod(tagger, computeis, vec, is, listed);
416: } else if (listed) *listed = PETSC_FALSE;
417: PetscFunctionReturn(PETSC_SUCCESS);
418: }
420: PetscErrorCode VecTaggerComputeIS_FromBoxes(VecTagger tagger, Vec vec, IS *is, PetscBool *listed)
421: {
422: PetscInt numBoxes;
423: VecTaggerBox *boxes;
424: PetscInt numTagged, offset;
425: PetscInt *tagged;
426: PetscInt bs, b, i, j, k, n;
427: PetscBool invert;
428: const PetscScalar *vecArray;
429: PetscBool boxlisted;
431: PetscFunctionBegin;
432: PetscCall(VecTaggerGetBlockSize(tagger, &bs));
433: PetscCall(VecTaggerComputeBoxes(tagger, vec, &numBoxes, &boxes, &boxlisted));
434: if (!boxlisted) {
435: if (listed) *listed = PETSC_FALSE;
436: PetscFunctionReturn(PETSC_SUCCESS);
437: }
438: PetscCall(VecGetArrayRead(vec, &vecArray));
439: PetscCall(VecGetLocalSize(vec, &n));
440: invert = tagger->invert;
441: numTagged = 0;
442: offset = 0;
443: tagged = NULL;
444: PetscCheck(n % bs == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "blocksize %" PetscInt_FMT " does not divide vector length %" PetscInt_FMT, bs, n);
445: n /= bs;
446: for (i = 0; i < 2; i++) {
447: if (i) PetscCall(PetscMalloc1(numTagged, &tagged));
448: for (j = 0; j < n; j++) {
449: for (k = 0; k < numBoxes; k++) {
450: for (b = 0; b < bs; b++) {
451: PetscScalar val = vecArray[j * bs + b];
452: PetscInt l = k * bs + b;
453: VecTaggerBox box;
454: PetscBool in;
456: box = boxes[l];
457: #if !defined(PETSC_USE_COMPLEX)
458: in = (PetscBool)((box.min <= val) && (val <= box.max));
459: #else
460: in = (PetscBool)((PetscRealPart(box.min) <= PetscRealPart(val)) && (PetscImaginaryPart(box.min) <= PetscImaginaryPart(val)) && (PetscRealPart(val) <= PetscRealPart(box.max)) && (PetscImaginaryPart(val) <= PetscImaginaryPart(box.max)));
461: #endif
462: if (!in) break;
463: }
464: if (b == bs) break;
465: }
466: if ((PetscBool)(k < numBoxes) ^ invert) {
467: if (!i) numTagged++;
468: else tagged[offset++] = j;
469: }
470: }
471: }
472: PetscCall(VecRestoreArrayRead(vec, &vecArray));
473: PetscCall(PetscFree(boxes));
474: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)vec), numTagged, tagged, PETSC_OWN_POINTER, is));
475: PetscCall(ISSort(*is));
476: if (listed) *listed = PETSC_TRUE;
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }