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: }