Actual source code: tagger.c

  1: #include <petsc/private/vecimpl.h>

  3: /*@C
  4:    VecTaggerCreate - create a `VecTagger` 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 compute 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 `VecTagger` will operate

 25:    Output Parameter:
 26: .  tagger - new Vec tagger context

 28:    Level: advanced

 30: .seealso: `VecTagger`, `VecTaggerSetBlockSize()`, `VecTaggerSetFromOptions()`, `VecTaggerSetUp()`, `VecTaggerComputeIS()`, `VecTaggerComputeBoxes()`, `VecTaggerDestroy()`
 31: @*/
 32: PetscErrorCode VecTaggerCreate(MPI_Comm comm, VecTagger *tagger)
 33: {
 34:   VecTagger b;

 36:   PetscFunctionBegin;
 38:   PetscCall(VecTaggerInitializePackage());

 40:   PetscCall(PetscHeaderCreate(b, VEC_TAGGER_CLASSID, "VecTagger", "Vec Tagger", "Vec", comm, VecTaggerDestroy, VecTaggerView));

 42:   b->blocksize   = 1;
 43:   b->invert      = PETSC_FALSE;
 44:   b->setupcalled = PETSC_FALSE;

 46:   *tagger = b;
 47:   PetscFunctionReturn(PETSC_SUCCESS);
 48: }

 50: /*@C
 51:    VecTaggerSetType - set the Vec tagger implementation

 53:    Collective

 55:    Input Parameters:
 56: +  tagger - the `VecTagger` context
 57: -  type - a known method

 59:    Options Database Key:
 60: .  -vec_tagger_type <type> - Sets the method; use -help for a list
 61:    of available methods (for instance, absolute, relative, cdf, or, and)

 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;

 84:   PetscCall(PetscObjectTypeCompare((PetscObject)tagger, type, &match));
 85:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

 87:   PetscCall(PetscFunctionListFind(VecTaggerList, type, &r));
 88:   PetscCheck(r, PETSC_COMM_SELF, 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: /*@C
 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;
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()`, `VecTaggerSetUp()`
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()`

188: @*/
189: PetscErrorCode VecTaggerSetFromOptions(VecTagger tagger)
190: {
191:   VecTaggerType deft;
192:   char          type[256];
193:   PetscBool     flg;

195:   PetscFunctionBegin;
197:   PetscObjectOptionsBegin((PetscObject)tagger);
198:   deft = ((PetscObject)tagger)->type_name ? ((PetscObject)tagger)->type_name : VECTAGGERABSOLUTE;
199:   PetscCall(PetscOptionsFList("-vec_tagger_type", "VecTagger implementation type", "VecTaggerSetType", VecTaggerList, deft, type, 256, &flg));
200:   PetscCall(VecTaggerSetType(tagger, flg ? type : deft));
201:   PetscCall(PetscOptionsInt("-vec_tagger_block_size", "block size of the vectors the tagger operates on", "VecTaggerSetBlockSize", tagger->blocksize, &tagger->blocksize, NULL));
202:   PetscCall(PetscOptionsBool("-vec_tagger_invert", "invert the set of indices returned by VecTaggerComputeIS()", "VecTaggerSetInvert", tagger->invert, &tagger->invert, NULL));
203:   PetscTryTypeMethod(tagger, setfromoptions, PetscOptionsObject);
204:   PetscOptionsEnd();
205:   PetscFunctionReturn(PETSC_SUCCESS);
206: }

208: /*@C
209:    VecTaggerSetBlockSize - block size of the set of indices returned by `VecTaggerComputeIS()`.  Values greater than one
210:    are useful when there are multiple criteria for determining which indices to include in the set.  For example,
211:    consider adaptive mesh refinement in a multiphysics problem, with metrics of solution quality for multiple fields
212:    measure on each cell.  The size of the vector will be [numCells * numFields]; the VecTagger block size should be
213:    numFields; VecTaggerComputeIS() will return indices in the range [0,numCells), i.e., one index is given for each
214:    block of values.

216:    Note that the block size of the vector does not have to match.

218:    Note also that the index set created in `VecTaggerComputeIS()` has block size: it is an index set over the list of
219:    items that the vector refers to, not to the vector itself.

221:    Logically Collective

223:    Input Parameters:
224: +  tagger - vec tagger
225: -  blocksize - block size of the criteria used to tagger vectors

227:    Level: advanced

229: .seealso: `VecTaggerComputeIS()`, `VecTaggerGetBlockSize()`, `VecSetBlockSize()`, `VecGetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
230: @*/
231: PetscErrorCode VecTaggerSetBlockSize(VecTagger tagger, PetscInt blocksize)
232: {
233:   PetscFunctionBegin;
236:   tagger->blocksize = blocksize;
237:   PetscFunctionReturn(PETSC_SUCCESS);
238: }

240: /*@C
241:    VecTaggerGetBlockSize - get the block size of the indices created by `VecTaggerComputeIS()`.

243:    Logically Collective

245:    Input Parameter:
246: .  tagger - vec tagger

248:    Output Parameter:
249: .  blocksize - block size of the vectors the tagger operates on

251:    Level: advanced

253: .seealso: `VecTaggerComputeIS()`, `VecTaggerSetBlockSize()`, `VecTagger`, `VecTaggerCreate()`
254: @*/
255: PetscErrorCode VecTaggerGetBlockSize(VecTagger tagger, PetscInt *blocksize)
256: {
257:   PetscFunctionBegin;
260:   *blocksize = tagger->blocksize;
261:   PetscFunctionReturn(PETSC_SUCCESS);
262: }

264: /*@C
265:    VecTaggerSetInvert - If the tagged index sets are based on boxes that can be returned by `VecTaggerComputeBoxes()`,
266:    then this option inverts values used to compute the IS, i.e., from being in the union of the boxes to being in the
267:    intersection of their exteriors.

269:    Logically Collective

271:    Input Parameters:
272: +  tagger - vec tagger
273: -  invert - `PETSC_TRUE` to invert, `PETSC_FALSE` to use the indices as is

275:    Level: advanced

277: .seealso: `VecTaggerComputeIS()`, `VecTaggerGetInvert()`, `VecTagger`, `VecTaggerCreate()`
278: @*/
279: PetscErrorCode VecTaggerSetInvert(VecTagger tagger, PetscBool invert)
280: {
281:   PetscFunctionBegin;
284:   tagger->invert = invert;
285:   PetscFunctionReturn(PETSC_SUCCESS);
286: }

288: /*@C
289:    VecTaggerGetInvert - get whether the set of indices returned by `VecTaggerComputeIS()` are inverted

291:    Logically Collective

293:    Input Parameter:
294: .  tagger - vec tagger

296:    Output Parameter:
297: .  invert - `PETSC_TRUE` to invert, `PETSC_FALSE` to use the indices as is

299:    Level: advanced

301: .seealso: `VecTaggerComputeIS()`, `VecTaggerSetInvert()`, `VecTagger`, `VecTaggerCreate()`
302: @*/
303: PetscErrorCode VecTaggerGetInvert(VecTagger tagger, PetscBool *invert)
304: {
305:   PetscFunctionBegin;
308:   *invert = tagger->invert;
309:   PetscFunctionReturn(PETSC_SUCCESS);
310: }

312: /*@C
313:    VecTaggerView - view a `VecTagger` context

315:    Collective

317:    Input Parameters:
318: +  tagger - vec tagger
319: -  viewer - viewer to display tagger, for example `PETSC_VIEWER_STDOUT_WORLD`

321:    Level: advanced

323: .seealso: `VecTaggerCreate()`, `VecTagger`, `VecTaggerCreate()`
324: @*/
325: PetscErrorCode VecTaggerView(VecTagger tagger, PetscViewer viewer)
326: {
327:   PetscBool iascii;

329:   PetscFunctionBegin;
331:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tagger), &viewer));
333:   PetscCheckSameComm(tagger, 1, viewer, 2);
334:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
335:   if (iascii) {
336:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)tagger, viewer));
337:     PetscCall(PetscViewerASCIIPushTab(viewer));
338:     PetscCall(PetscViewerASCIIPrintf(viewer, "Block size: %" PetscInt_FMT "\n", tagger->blocksize));
339:     PetscTryTypeMethod(tagger, view, viewer);
340:     if (tagger->invert) PetscCall(PetscViewerASCIIPrintf(viewer, "Inverting ISs.\n"));
341:     PetscCall(PetscViewerASCIIPopTab(viewer));
342:   }
343:   PetscFunctionReturn(PETSC_SUCCESS);
344: }

346: /*@C
347:    VecTaggerComputeBoxes - If the tagged index set can be summarized as a list of boxes of values, returns that list, otherwise returns
348:          in listed `PETSC_FALSE`

350:    Collective on VecTagger

352:    Input Parameters:
353: +  tagger - the `VecTagger` context
354: -  vec - the vec to tag

356:    Output Parameters:
357: +  numBoxes - the number of boxes in the tag definition
358: .  boxes - a newly allocated list of boxes.  This is a flat array of (BlockSize * `numBoxe`s) pairs that the user can free with `PetscFree()`.
359: -  listed - `PETSC_TRUE` if a list was created, pass in `NULL` if not needed

361:    Level: advanced

363:    Note:
364:      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.

366: .seealso: `VecTaggerComputeIS()`, `VecTagger`, `VecTaggerCreate()`
367: @*/
368: PetscErrorCode VecTaggerComputeBoxes(VecTagger tagger, Vec vec, PetscInt *numBoxes, VecTaggerBox **boxes, PetscBool *listed)
369: {
370:   PetscInt vls, tbs;

372:   PetscFunctionBegin;
377:   PetscCall(VecGetLocalSize(vec, &vls));
378:   PetscCall(VecTaggerGetBlockSize(tagger, &tbs));
379:   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);
380:   if (tagger->ops->computeboxes) {
381:     *listed = PETSC_TRUE;
382:     PetscCall((*tagger->ops->computeboxes)(tagger, vec, numBoxes, boxes, listed));
383:   } else *listed = PETSC_FALSE;
384:   PetscFunctionReturn(PETSC_SUCCESS);
385: }

387: /*@C
388:    VecTaggerComputeIS - Use a `VecTagger` context to tag a set of indices based on a vector's values

390:    Collective

392:    Input Parameters:
393: +  tagger - the `VecTagger` context
394: -  vec - the vec to tag

396:    Output Parameters:
397: +  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()`.
398: -  listed - routine was able to compute the `IS`, pass in `NULL` if not needed

400:    Level: advanced

402: .seealso: `VecTaggerComputeBoxes()`, `VecTagger`, `VecTaggerCreate()`
403: @*/
404: PetscErrorCode VecTaggerComputeIS(VecTagger tagger, Vec vec, IS *is, PetscBool *listed)
405: {
406:   PetscInt vls, tbs;

408:   PetscFunctionBegin;
412:   PetscCall(VecGetLocalSize(vec, &vls));
413:   PetscCall(VecTaggerGetBlockSize(tagger, &tbs));
414:   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);
415:   if (tagger->ops->computeis) {
416:     PetscCall((*tagger->ops->computeis)(tagger, vec, is, listed));
417:   } else if (listed) *listed = PETSC_FALSE;
418:   PetscFunctionReturn(PETSC_SUCCESS);
419: }

421: PetscErrorCode VecTaggerComputeIS_FromBoxes(VecTagger tagger, Vec vec, IS *is, PetscBool *listed)
422: {
423:   PetscInt           numBoxes;
424:   VecTaggerBox      *boxes;
425:   PetscInt           numTagged, offset;
426:   PetscInt          *tagged;
427:   PetscInt           bs, b, i, j, k, n;
428:   PetscBool          invert;
429:   const PetscScalar *vecArray;
430:   PetscBool          boxlisted;

432:   PetscFunctionBegin;
433:   PetscCall(VecTaggerGetBlockSize(tagger, &bs));
434:   PetscCall(VecTaggerComputeBoxes(tagger, vec, &numBoxes, &boxes, &boxlisted));
435:   if (!boxlisted) {
436:     if (listed) *listed = PETSC_FALSE;
437:     PetscFunctionReturn(PETSC_SUCCESS);
438:   }
439:   PetscCall(VecGetArrayRead(vec, &vecArray));
440:   PetscCall(VecGetLocalSize(vec, &n));
441:   invert    = tagger->invert;
442:   numTagged = 0;
443:   offset    = 0;
444:   tagged    = NULL;
445:   PetscCheck(n % bs == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "blocksize %" PetscInt_FMT " does not divide vector length %" PetscInt_FMT, bs, n);
446:   n /= bs;
447:   for (i = 0; i < 2; i++) {
448:     if (i) PetscCall(PetscMalloc1(numTagged, &tagged));
449:     for (j = 0; j < n; j++) {
450:       for (k = 0; k < numBoxes; k++) {
451:         for (b = 0; b < bs; b++) {
452:           PetscScalar  val = vecArray[j * bs + b];
453:           PetscInt     l   = k * bs + b;
454:           VecTaggerBox box;
455:           PetscBool    in;

457:           box = boxes[l];
458: #if !defined(PETSC_USE_COMPLEX)
459:           in = (PetscBool)((box.min <= val) && (val <= box.max));
460: #else
461:           in = (PetscBool)((PetscRealPart(box.min) <= PetscRealPart(val)) && (PetscImaginaryPart(box.min) <= PetscImaginaryPart(val)) && (PetscRealPart(val) <= PetscRealPart(box.max)) && (PetscImaginaryPart(val) <= PetscImaginaryPart(box.max)));
462: #endif
463:           if (!in) break;
464:         }
465:         if (b == bs) break;
466:       }
467:       if ((PetscBool)(k < numBoxes) ^ invert) {
468:         if (!i) numTagged++;
469:         else tagged[offset++] = j;
470:       }
471:     }
472:   }
473:   PetscCall(VecRestoreArrayRead(vec, &vecArray));
474:   PetscCall(PetscFree(boxes));
475:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)vec), numTagged, tagged, PETSC_OWN_POINTER, is));
476:   PetscCall(ISSort(*is));
477:   if (listed) *listed = PETSC_TRUE;
478:   PetscFunctionReturn(PETSC_SUCCESS);
479: }