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:   if (tagger->ops->destroy) (*(tagger)->ops->destroy)(tagger);
 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) {*tagger = NULL; return 0;}
137:   if ((*tagger)->ops->destroy) (*(*tagger)->ops->destroy)(*tagger);
138:   PetscHeaderDestroy(tagger);
139:   return 0;
140: }

142: /*@
143:    VecTaggerSetUp - set up a VecTagger context

145:    Collective

147:    Input Parameter:
148: .  tagger - Vec tagger object

150:    Level: advanced

152: .seealso: VecTaggerSetFromOptions(), VecTaggerSetType(), VecTagger, VecTaggerCreate(), VecTaggerSetUp()
153: @*/
154: PetscErrorCode VecTaggerSetUp(VecTagger tagger)
155: {
156:   if (tagger->setupcalled) return 0;
157:   if (!((PetscObject)tagger)->type_name) VecTaggerSetType(tagger,VECTAGGERABSOLUTE);
158:   if (tagger->ops->setup) (*tagger->ops->setup)(tagger);
159:   tagger->setupcalled = PETSC_TRUE;
160:   return 0;
161: }

163: /*@C
164:    VecTaggerSetFromOptions - set VecTagger options using the options database

166:    Logically Collective

168:    Input Parameter:
169: .  tagger - vec tagger

171:    Options Database Keys:
172: +  -vec_tagger_type       - implementation type, see VecTaggerSetType()
173: .  -vec_tagger_block_size - set the block size, see VecTaggerSetBlockSize()
174: -  -vec_tagger_invert     - invert the index set returned by VecTaggerComputeIS()

176:    Level: advanced

178: .seealso: VecTagger, VecTaggerCreate(), VecTaggerSetUp()

180: @*/
181: PetscErrorCode VecTaggerSetFromOptions(VecTagger tagger)
182: {
183:   VecTaggerType  deft;
184:   char           type[256];
186:   PetscBool      flg;

189:   PetscObjectOptionsBegin((PetscObject)tagger);
190:   deft = ((PetscObject)tagger)->type_name ? ((PetscObject)tagger)->type_name : VECTAGGERABSOLUTE;
191:   PetscOptionsFList("-vec_tagger_type","VecTagger implementation type","VecTaggerSetType",VecTaggerList,deft,type,256,&flg);
192:   VecTaggerSetType(tagger,flg ? type : deft);
193:   PetscOptionsInt("-vec_tagger_block_size","block size of the vectors the tagger operates on","VecTaggerSetBlockSize",tagger->blocksize,&tagger->blocksize,NULL);
194:   PetscOptionsBool("-vec_tagger_invert","invert the set of indices returned by VecTaggerComputeIS()","VecTaggerSetInvert",tagger->invert,&tagger->invert,NULL);
195:   if (tagger->ops->setfromoptions) (*tagger->ops->setfromoptions)(PetscOptionsObject,tagger);
196:   PetscOptionsEnd();
197:   return 0;
198: }

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

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

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

213:    Logically Collective

215:    Input Parameters:
216: +  tagger - vec tagger
217: -  blocksize - block size of the criteria used to tagger vectors

219:    Level: advanced

221: .seealso: VecTaggerComputeIS(), VecTaggerGetBlockSize(), VecSetBlockSize(), VecGetBlockSize(), VecTagger, VecTaggerCreate()
222: @*/
223: PetscErrorCode VecTaggerSetBlockSize(VecTagger tagger, PetscInt blocksize)
224: {
227:   tagger->blocksize = blocksize;
228:   return 0;
229: }

231: /*@C
232:    VecTaggerGetBlockSize - get the block size of the indices created by VecTaggerComputeIS().

234:    Logically Collective

236:    Input Parameter:
237: .  tagger - vec tagger

239:    Output Parameter:
240: .  blocksize - block size of the vectors the tagger operates on

242:    Level: advanced

244: .seealso: VecTaggerComputeIS(), VecTaggerSetBlockSize(), VecTagger, VecTaggerCreate()
245: @*/
246: PetscErrorCode VecTaggerGetBlockSize(VecTagger tagger, PetscInt *blocksize)
247: {
250:   *blocksize = tagger->blocksize;
251:   return 0;
252: }

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

259:    Logically Collective

261:    Input Parameters:
262: +  tagger - vec tagger
263: -  invert - PETSC_TRUE to invert, PETSC_FALSE to use the indices as is

265:    Level: advanced

267: .seealso: VecTaggerComputeIS(), VecTaggerGetInvert(), VecTagger, VecTaggerCreate()
268: @*/
269: PetscErrorCode VecTaggerSetInvert(VecTagger tagger, PetscBool invert)
270: {
273:   tagger->invert = invert;
274:   return 0;
275: }

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

280:    Logically Collective

282:    Input Parameter:
283: .  tagger - vec tagger

285:    Output Parameter:
286: .  invert - PETSC_TRUE to invert, PETSC_FALSE to use the indices as is

288:    Level: advanced

290: .seealso: VecTaggerComputeIS(), VecTaggerSetInvert(), VecTagger, VecTaggerCreate()
291: @*/
292: PetscErrorCode VecTaggerGetInvert(VecTagger tagger, PetscBool *invert)
293: {
296:   *invert = tagger->invert;
297:   return 0;
298: }

300: /*@C
301:    VecTaggerView - view a VecTagger context

303:    Collective

305:    Input Parameters:
306: +  tagger - vec tagger
307: -  viewer - viewer to display tagger, for example PETSC_VIEWER_STDOUT_WORLD

309:    Level: advanced

311: .seealso: VecTaggerCreate(), VecTagger, VecTaggerCreate()
312: @*/
313: PetscErrorCode VecTaggerView(VecTagger tagger,PetscViewer viewer)
314: {
315:   PetscBool         iascii;

318:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tagger),&viewer);
321:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
322:   if (iascii) {
323:     PetscObjectPrintClassNamePrefixType((PetscObject)tagger,viewer);
324:     PetscViewerASCIIPushTab(viewer);
325:     PetscViewerASCIIPrintf(viewer,"Block size: %" PetscInt_FMT "\n",tagger->blocksize);
326:     if (tagger->ops->view) (*tagger->ops->view)(tagger,viewer);
327:     if (tagger->invert) PetscViewerASCIIPrintf(viewer,"Inverting ISs.\n");
328:     PetscViewerASCIIPopTab(viewer);
329:   }
330:   return 0;
331: }

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

337:    Collective on VecTagger

339:    Input Parameters:
340: +  tagger - the VecTagger context
341: -  vec - the vec to tag

343:    Output Parameters:
344: +  numBoxes - the number of boxes in the tag definition
345: .  boxes - a newly allocated list of boxes.  This is a flat array of (BlockSize * numBoxes) pairs that the user can free with PetscFree().
346: -  listed - PETSC_TRUE if a list was created, pass in NULL if not needed

348:    Notes:
349:      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.

351:    Level: advanced

353: .seealso: VecTaggerComputeIS(), VecTagger, VecTaggerCreate()
354: @*/
355: PetscErrorCode VecTaggerComputeBoxes(VecTagger tagger,Vec vec,PetscInt *numBoxes,VecTaggerBox **boxes,PetscBool *listed)
356: {
357:   PetscInt       vls, tbs;

363:   VecGetLocalSize(vec,&vls);
364:   VecTaggerGetBlockSize(tagger,&tbs);
366:   if (tagger->ops->computeboxes) {
367:     *listed = PETSC_TRUE;
368:     (*tagger->ops->computeboxes) (tagger,vec,numBoxes,boxes,listed);
369:   }  else *listed = PETSC_FALSE;
370:   return 0;
371: }

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

376:    Collective on VecTagger

378:    Input Parameters:
379: +  tagger - the VecTagger context
380: -  vec - the vec to tag

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

386:    Level: advanced

388: .seealso: VecTaggerComputeBoxes(), VecTagger, VecTaggerCreate()
389: @*/
390: PetscErrorCode VecTaggerComputeIS(VecTagger tagger,Vec vec,IS *is,PetscBool *listed)
391: {
392:   PetscInt       vls, tbs;

397:   VecGetLocalSize(vec,&vls);
398:   VecTaggerGetBlockSize(tagger,&tbs);
400:   if (tagger->ops->computeis) {
401:     (*tagger->ops->computeis) (tagger,vec,is,listed);
402:   } else if (listed) *listed = PETSC_FALSE;
403:   return 0;
404: }

406: PetscErrorCode VecTaggerComputeIS_FromBoxes(VecTagger tagger, Vec vec, IS *is,PetscBool *listed)
407: {
408:   PetscInt          numBoxes;
409:   VecTaggerBox      *boxes;
410:   PetscInt          numTagged, offset;
411:   PetscInt          *tagged;
412:   PetscInt          bs, b, i, j, k, n;
413:   PetscBool         invert;
414:   const PetscScalar *vecArray;
415:   PetscBool         boxlisted;

417:   VecTaggerGetBlockSize(tagger,&bs);
418:   VecTaggerComputeBoxes(tagger,vec,&numBoxes,&boxes,&boxlisted);
419:   if (!boxlisted) {
420:     if (listed) *listed = PETSC_FALSE;
421:     return 0;
422:   }
423:   VecGetArrayRead(vec, &vecArray);
424:   VecGetLocalSize(vec, &n);
425:   invert = tagger->invert;
426:   numTagged = 0;
427:   offset = 0;
428:   tagged = NULL;
430:   n /= bs;
431:   for (i = 0; i < 2; i++) {
432:     if (i) {
433:       PetscMalloc1(numTagged,&tagged);
434:     }
435:     for (j = 0; j < n; j++) {

437:       for (k = 0; k < numBoxes; k++) {
438:         for (b = 0; b < bs; b++) {
439:           PetscScalar  val = vecArray[j * bs + b];
440:           PetscInt     l = k * bs + b;
441:           VecTaggerBox box;
442:           PetscBool    in;

444:           box = boxes[l];
445: #if !defined(PETSC_USE_COMPLEX)
446:           in = (PetscBool) ((box.min <= val) && (val <= box.max));
447: #else
448:           in = (PetscBool) ((PetscRealPart     (box.min) <= PetscRealPart     (val)) &&
449:                             (PetscImaginaryPart(box.min) <= PetscImaginaryPart(val)) &&
450:                             (PetscRealPart     (val)     <= PetscRealPart     (box.max)) &&
451:                             (PetscImaginaryPart(val)     <= PetscImaginaryPart(box.max)));
452: #endif
453:           if (!in) break;
454:         }
455:         if (b == bs) break;
456:       }
457:       if ((PetscBool)(k < numBoxes) ^ invert) {
458:         if (!i) numTagged++;
459:         else    tagged[offset++] = j;
460:       }
461:     }
462:   }
463:   VecRestoreArrayRead(vec, &vecArray);
464:   PetscFree(boxes);
465:   ISCreateGeneral(PetscObjectComm((PetscObject)vec),numTagged,tagged,PETSC_OWN_POINTER,is);
466:   ISSort(*is);
467:   if (listed) *listed = PETSC_TRUE;
468:   return 0;
469: }