Actual source code: tagger.c
1: #include <petsc/private/vecimpl.h>
3: /*@C
4: VecTaggerCreate - create a Vec tagger context. This object is used to control the tagging/selection of index sets
5: based on the values in a vector. This is used, for example, in adaptive simulations when aspects are selected for
6: refinement or coarsening. The primary intent is that the selected index sets are based purely on the values in the
7: vector, though implementations that do not follow this intent are possible.
9: Once a VecTagger is created (VecTaggerCreate()), optionally modified by options (VecTaggerSetFromOptions()), and
10: set up (VecTaggerSetUp()), it is applied to vectors with VecTaggerComputeIS() to comute the selected index sets.
12: In many cases, the selection criteria for an index is whether the corresponding value falls within a collection of
13: boxes: for this common case, VecTaggerCreateBoxes() can also be used to determine those boxes.
15: Provided implementations support tagging based on a box/interval of values (VECTAGGERABSOLUTE), based on a box of
16: values of relative to the range of values present in the vector (VECTAGGERRELATIVE), based on where values fall in
17: the cumulative distribution of values in the vector (VECTAGGERCDF), and based on unions (VECTAGGEROR) or
18: intersections (VECTAGGERAND) of other criteria.
20: Collective
22: Input Parameter:
23: . comm - communicator on which the vec tagger will operate
25: Output Parameter:
26: . tagger - new Vec tagger context
28: Level: advanced
30: .seealso: `VecTaggerSetBlockSize()`, `VecTaggerSetFromOptions()`, `VecTaggerSetUp()`, `VecTaggerComputeIS()`, `VecTaggerComputeBoxes()`, `VecTaggerDestroy()`
31: @*/
32: PetscErrorCode VecTaggerCreate(MPI_Comm comm, VecTagger *tagger)
33: {
34: VecTagger b;
37: VecTaggerInitializePackage();
39: PetscHeaderCreate(b, VEC_TAGGER_CLASSID, "VecTagger", "Vec Tagger", "Vec", comm, VecTaggerDestroy, VecTaggerView);
41: b->blocksize = 1;
42: b->invert = PETSC_FALSE;
43: b->setupcalled = PETSC_FALSE;
45: *tagger = b;
46: return 0;
47: }
49: /*@C
50: VecTaggerSetType - set the Vec tagger implementation
52: Collective on VecTagger
54: Input Parameters:
55: + tagger - the VecTagger context
56: - type - a known method
58: Options Database Key:
59: . -vec_tagger_type <type> - Sets the method; use -help for a list
60: of available methods (for instance, absolute, relative, cdf, or, and)
62: Notes:
63: See "include/petscvec.h" for available methods (for instance)
64: + VECTAGGERABSOLUTE - tag based on a box of values
65: . VECTAGGERRELATIVE - tag based on a box relative to the range of values present in the vector
66: . VECTAGGERCDF - tag based on a box in the cumulative distribution of values present in the vector
67: . VECTAGGEROR - tag based on the union of a set of VecTagger contexts
68: - VECTAGGERAND - tag based on the intersection of a set of other VecTagger contexts
70: Level: advanced
72: .seealso: `VecTaggerType`, `VecTaggerCreate()`, `VecTagger`
73: @*/
74: PetscErrorCode VecTaggerSetType(VecTagger tagger, VecTaggerType type)
75: {
76: PetscBool match;
77: PetscErrorCode (*r)(VecTagger);
82: PetscObjectTypeCompare((PetscObject)tagger, type, &match);
83: if (match) return 0;
85: PetscFunctionListFind(VecTaggerList, type, &r);
87: /* Destroy the previous private VecTagger context */
88: PetscTryTypeMethod(tagger, destroy);
89: PetscMemzero(tagger->ops, sizeof(*tagger->ops));
90: PetscObjectChangeTypeName((PetscObject)tagger, type);
91: tagger->ops->create = r;
92: (*r)(tagger);
93: return 0;
94: }
96: /*@C
97: VecTaggerGetType - Gets the VecTagger type name (as a string) from the VecTagger.
99: Not Collective
101: Input Parameter:
102: . tagger - The Vec tagger context
104: Output Parameter:
105: . type - The VecTagger type name
107: Level: advanced
109: .seealso: `VecTaggerSetType()`, `VecTaggerCreate()`, `VecTaggerSetFromOptions()`, `VecTagger`, `VecTaggerType`
110: @*/
111: PetscErrorCode VecTaggerGetType(VecTagger tagger, VecTaggerType *type)
112: {
115: VecTaggerRegisterAll();
116: *type = ((PetscObject)tagger)->type_name;
117: return 0;
118: }
120: /*@
121: VecTaggerDestroy - destroy a VecTagger context
123: Collective
125: Input Parameter:
126: . tagger - address of tagger
128: Level: advanced
130: .seealso: `VecTaggerCreate()`, `VecTaggerSetType()`, `VecTagger`
131: @*/
132: PetscErrorCode VecTaggerDestroy(VecTagger *tagger)
133: {
134: if (!*tagger) return 0;
136: if (--((PetscObject)(*tagger))->refct > 0) {
137: *tagger = NULL;
138: return 0;
139: }
140: PetscTryTypeMethod((*tagger), destroy);
141: PetscHeaderDestroy(tagger);
142: return 0;
143: }
145: /*@
146: VecTaggerSetUp - set up a VecTagger context
148: Collective
150: Input Parameter:
151: . tagger - Vec tagger object
153: Level: advanced
155: .seealso: `VecTaggerSetFromOptions()`, `VecTaggerSetType()`, `VecTagger`, `VecTaggerCreate()`, `VecTaggerSetUp()`
156: @*/
157: PetscErrorCode VecTaggerSetUp(VecTagger tagger)
158: {
159: if (tagger->setupcalled) return 0;
160: if (!((PetscObject)tagger)->type_name) VecTaggerSetType(tagger, VECTAGGERABSOLUTE);
161: PetscTryTypeMethod(tagger, setup);
162: tagger->setupcalled = PETSC_TRUE;
163: return 0;
164: }
166: /*@C
167: VecTaggerSetFromOptions - set VecTagger options using the options database
169: Logically Collective
171: Input Parameter:
172: . tagger - vec tagger
174: Options Database Keys:
175: + -vec_tagger_type - implementation type, see VecTaggerSetType()
176: . -vec_tagger_block_size - set the block size, see VecTaggerSetBlockSize()
177: - -vec_tagger_invert - invert the index set returned by VecTaggerComputeIS()
179: Level: advanced
181: .seealso: `VecTagger`, `VecTaggerCreate()`, `VecTaggerSetUp()`
183: @*/
184: PetscErrorCode VecTaggerSetFromOptions(VecTagger tagger)
185: {
186: VecTaggerType deft;
187: char type[256];
188: PetscBool flg;
191: PetscObjectOptionsBegin((PetscObject)tagger);
192: deft = ((PetscObject)tagger)->type_name ? ((PetscObject)tagger)->type_name : VECTAGGERABSOLUTE;
193: PetscOptionsFList("-vec_tagger_type", "VecTagger implementation type", "VecTaggerSetType", VecTaggerList, deft, type, 256, &flg);
194: VecTaggerSetType(tagger, flg ? type : deft);
195: PetscOptionsInt("-vec_tagger_block_size", "block size of the vectors the tagger operates on", "VecTaggerSetBlockSize", tagger->blocksize, &tagger->blocksize, NULL);
196: PetscOptionsBool("-vec_tagger_invert", "invert the set of indices returned by VecTaggerComputeIS()", "VecTaggerSetInvert", tagger->invert, &tagger->invert, NULL);
197: PetscTryTypeMethod(tagger, setfromoptions, PetscOptionsObject);
198: PetscOptionsEnd();
199: return 0;
200: }
202: /*@C
203: VecTaggerSetBlockSize - block size of the set of indices returned by VecTaggerComputeIS(). Values greater than one
204: are useful when there are multiple criteria for determining which indices to include in the set. For example,
205: consider adaptive mesh refinement in a multiphysics problem, with metrics of solution quality for multiple fields
206: measure on each cell. The size of the vector will be [numCells * numFields]; the VecTagger block size should be
207: numFields; VecTaggerComputeIS() will return indices in the range [0,numCells), i.e., one index is given for each
208: block of values.
210: Note that the block size of the vector does not have to match.
212: Note also that the index set created in VecTaggerComputeIS() has block size: it is an index set over the list of
213: items that the vector refers to, not to the vector itself.
215: Logically Collective
217: Input Parameters:
218: + tagger - vec tagger
219: - blocksize - block size of the criteria used to tagger vectors
221: Level: advanced
223: .seealso: `VecTaggerComputeIS()`, `VecTaggerGetBlockSize()`, `VecSetBlockSize()`, `VecGetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
224: @*/
225: PetscErrorCode VecTaggerSetBlockSize(VecTagger tagger, PetscInt blocksize)
226: {
229: tagger->blocksize = blocksize;
230: return 0;
231: }
233: /*@C
234: VecTaggerGetBlockSize - get the block size of the indices created by VecTaggerComputeIS().
236: Logically Collective
238: Input Parameter:
239: . tagger - vec tagger
241: Output Parameter:
242: . blocksize - block size of the vectors the tagger operates on
244: Level: advanced
246: .seealso: `VecTaggerComputeIS()`, `VecTaggerSetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
247: @*/
248: PetscErrorCode VecTaggerGetBlockSize(VecTagger tagger, PetscInt *blocksize)
249: {
252: *blocksize = tagger->blocksize;
253: return 0;
254: }
256: /*@C
257: VecTaggerSetInvert - If the tagged index sets are based on boxes that can be returned by VecTaggerComputeBoxes(),
258: then this option inverts values used to compute the IS, i.e., from being in the union of the boxes to being in the
259: intersection of their exteriors.
261: Logically Collective
263: Input Parameters:
264: + tagger - vec tagger
265: - invert - PETSC_TRUE to invert, PETSC_FALSE to use the indices as is
267: Level: advanced
269: .seealso: `VecTaggerComputeIS()`, `VecTaggerGetInvert()`, `VecTagger`, `VecTaggerCreate()`
270: @*/
271: PetscErrorCode VecTaggerSetInvert(VecTagger tagger, PetscBool invert)
272: {
275: tagger->invert = invert;
276: return 0;
277: }
279: /*@C
280: VecTaggerGetInvert - get whether the set of indices returned by VecTaggerComputeIS() are inverted
282: Logically Collective
284: Input Parameter:
285: . tagger - vec tagger
287: Output Parameter:
288: . invert - PETSC_TRUE to invert, PETSC_FALSE to use the indices as is
290: Level: advanced
292: .seealso: `VecTaggerComputeIS()`, `VecTaggerSetInvert()`, `VecTagger`, `VecTaggerCreate()`
293: @*/
294: PetscErrorCode VecTaggerGetInvert(VecTagger tagger, PetscBool *invert)
295: {
298: *invert = tagger->invert;
299: return 0;
300: }
302: /*@C
303: VecTaggerView - view a VecTagger context
305: Collective
307: Input Parameters:
308: + tagger - vec tagger
309: - viewer - viewer to display tagger, for example PETSC_VIEWER_STDOUT_WORLD
311: Level: advanced
313: .seealso: `VecTaggerCreate()`, `VecTagger`, `VecTaggerCreate()`
314: @*/
315: PetscErrorCode VecTaggerView(VecTagger tagger, PetscViewer viewer)
316: {
317: PetscBool iascii;
320: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tagger), &viewer);
323: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
324: if (iascii) {
325: PetscObjectPrintClassNamePrefixType((PetscObject)tagger, viewer);
326: PetscViewerASCIIPushTab(viewer);
327: PetscViewerASCIIPrintf(viewer, "Block size: %" PetscInt_FMT "\n", tagger->blocksize);
328: PetscTryTypeMethod(tagger, view, viewer);
329: if (tagger->invert) PetscViewerASCIIPrintf(viewer, "Inverting ISs.\n");
330: PetscViewerASCIIPopTab(viewer);
331: }
332: return 0;
333: }
335: /*@C
336: VecTaggerComputeBoxes - If the tagged index set can be summarized as a list of boxes of values, returns that list, otherwise returns
337: in listed PETSC_FALSE
339: Collective on VecTagger
341: Input Parameters:
342: + tagger - the VecTagger context
343: - vec - the vec to tag
345: Output Parameters:
346: + numBoxes - the number of boxes in the tag definition
347: . boxes - a newly allocated list of boxes. This is a flat array of (BlockSize * numBoxes) pairs that the user can free with PetscFree().
348: - listed - PETSC_TRUE if a list was created, pass in NULL if not needed
350: Notes:
351: 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.
353: Level: advanced
355: .seealso: `VecTaggerComputeIS()`, `VecTagger`, `VecTaggerCreate()`
356: @*/
357: PetscErrorCode VecTaggerComputeBoxes(VecTagger tagger, Vec vec, PetscInt *numBoxes, VecTaggerBox **boxes, PetscBool *listed)
358: {
359: PetscInt vls, tbs;
365: VecGetLocalSize(vec, &vls);
366: VecTaggerGetBlockSize(tagger, &tbs);
368: if (tagger->ops->computeboxes) {
369: *listed = PETSC_TRUE;
370: (*tagger->ops->computeboxes)(tagger, vec, numBoxes, boxes, listed);
371: } else *listed = PETSC_FALSE;
372: return 0;
373: }
375: /*@C
376: VecTaggerComputeIS - Use a VecTagger context to tag a set of indices based on a vector's values
378: Collective on VecTagger
380: Input Parameters:
381: + tagger - the VecTagger context
382: - vec - the vec to tag
384: Output Parameters:
385: + 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().
386: - listed - routine was able to compute the IS, pass in NULL if not needed
388: Level: advanced
390: .seealso: `VecTaggerComputeBoxes()`, `VecTagger`, `VecTaggerCreate()`
391: @*/
392: PetscErrorCode VecTaggerComputeIS(VecTagger tagger, Vec vec, IS *is, PetscBool *listed)
393: {
394: PetscInt vls, tbs;
399: VecGetLocalSize(vec, &vls);
400: VecTaggerGetBlockSize(tagger, &tbs);
402: if (tagger->ops->computeis) {
403: (*tagger->ops->computeis)(tagger, vec, is, listed);
404: } else if (listed) *listed = PETSC_FALSE;
405: return 0;
406: }
408: PetscErrorCode VecTaggerComputeIS_FromBoxes(VecTagger tagger, Vec vec, IS *is, PetscBool *listed)
409: {
410: PetscInt numBoxes;
411: VecTaggerBox *boxes;
412: PetscInt numTagged, offset;
413: PetscInt *tagged;
414: PetscInt bs, b, i, j, k, n;
415: PetscBool invert;
416: const PetscScalar *vecArray;
417: PetscBool boxlisted;
419: VecTaggerGetBlockSize(tagger, &bs);
420: VecTaggerComputeBoxes(tagger, vec, &numBoxes, &boxes, &boxlisted);
421: if (!boxlisted) {
422: if (listed) *listed = PETSC_FALSE;
423: return 0;
424: }
425: VecGetArrayRead(vec, &vecArray);
426: VecGetLocalSize(vec, &n);
427: invert = tagger->invert;
428: numTagged = 0;
429: offset = 0;
430: tagged = NULL;
432: n /= bs;
433: for (i = 0; i < 2; i++) {
434: if (i) PetscMalloc1(numTagged, &tagged);
435: for (j = 0; j < n; j++) {
436: for (k = 0; k < numBoxes; k++) {
437: for (b = 0; b < bs; b++) {
438: PetscScalar val = vecArray[j * bs + b];
439: PetscInt l = k * bs + b;
440: VecTaggerBox box;
441: PetscBool in;
443: box = boxes[l];
444: #if !defined(PETSC_USE_COMPLEX)
445: in = (PetscBool)((box.min <= val) && (val <= box.max));
446: #else
447: in = (PetscBool)((PetscRealPart(box.min) <= PetscRealPart(val)) && (PetscImaginaryPart(box.min) <= PetscImaginaryPart(val)) && (PetscRealPart(val) <= PetscRealPart(box.max)) && (PetscImaginaryPart(val) <= PetscImaginaryPart(box.max)));
448: #endif
449: if (!in) break;
450: }
451: if (b == bs) break;
452: }
453: if ((PetscBool)(k < numBoxes) ^ invert) {
454: if (!i) numTagged++;
455: else tagged[offset++] = j;
456: }
457: }
458: }
459: VecRestoreArrayRead(vec, &vecArray);
460: PetscFree(boxes);
461: ISCreateGeneral(PetscObjectComm((PetscObject)vec), numTagged, tagged, PETSC_OWN_POINTER, is);
462: ISSort(*is);
463: if (listed) *listed = PETSC_TRUE;
464: return 0;
465: }