Actual source code: tagger.c

petsc-master 2019-09-22
Report Typos and Errors
  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 Arguments:
 23: .  comm - communicator on which the vec tagger will operate

 25:    Output Arguments:
 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: {
 35:   VecTagger      b;

 39:   VecTaggerInitializePackage();

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

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

 47:   *tagger = b;
 48:   return(0);
 49: }

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

 54:    Collective on VecTagger

 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; use -help for a list
 62:    of available methods (for instance, absolute, relative, cdf, or, and)

 64:    Notes:
 65:    See "include/petscvec.h" for available methods (for instance)
 66: +    VECTAGGERABSOLUTE - tag based on a box of values
 67: .    VECTAGGERRELATIVE - tag based on a box relative to the range of values present in the vector
 68: .    VECTAGGERCDF      - tag based on a box in the cumulative distribution of values present in the vector
 69: .    VECTAGGEROR       - tag based on the union of a set of VecTagger contexts
 70: -    VECTAGGERAND      - tag based on the intersection of a set of other VecTagger contexts

 72:   Level: advanced

 74: .seealso: VecTaggerType, VecTaggerCreate()
 75: @*/
 76: PetscErrorCode VecTaggerSetType(VecTagger tagger,VecTaggerType type)
 77: {
 78:   PetscErrorCode ierr,(*r)(VecTagger);
 79:   PetscBool      match;


 85:   PetscObjectTypeCompare((PetscObject)tagger,type,&match);
 86:   if (match) return(0);

 88:   PetscFunctionListFind(VecTaggerList,type,&r);
 89:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested VecTagger type %s",type);
 90:   /* Destroy the previous private VecTagger context */
 91:   if (tagger->ops->destroy) {
 92:     (*(tagger)->ops->destroy)(tagger);
 93:   }
 94:   PetscMemzero(tagger->ops,sizeof(*tagger->ops));
 95:   PetscObjectChangeTypeName((PetscObject)tagger,type);
 96:   tagger->ops->create = r;
 97:   (*r)(tagger);
 98:   return(0);
 99: }

101: /*@C
102:   VecTaggerGetType - Gets the VecTagger type name (as a string) from the VecTagger.

104:   Not Collective

106:   Input Parameter:
107: . tagger  - The Vec tagger context

109:   Output Parameter:
110: . type - The VecTagger type name

112:   Level: advanced

114: .seealso: VecTaggerSetType(), VecTaggerCreate()
115: @*/
116: PetscErrorCode  VecTaggerGetType(VecTagger tagger, VecTaggerType *type)
117: {

123:   VecTaggerRegisterAll();
124:   *type = ((PetscObject)tagger)->type_name;
125:   return(0);
126: }

128: /*@
129:    VecTaggerDestroy - destroy a VecTagger context

131:    Collective

133:    Input Arguments:
134: .  tagger - address of tagger

136:    Level: advanced

138: .seealso: VecTaggerCreate()
139: @*/
140: PetscErrorCode VecTaggerDestroy(VecTagger *tagger)
141: {

145:   if (!*tagger) return(0);
147:   if (--((PetscObject)(*tagger))->refct > 0) {*tagger = 0; return(0);}
148:   if ((*tagger)->ops->destroy) {(*(*tagger)->ops->destroy)(*tagger);}
149:   PetscHeaderDestroy(tagger);
150:   return(0);
151: }

153: /*@
154:    VecTaggerSetUp - set up a VecTagger context

156:    Collective

158:    Input Arguments:
159: .  tagger - Vec tagger object

161:    Level: advanced

163: .seealso: VecTaggerSetFromOptions(), VecTaggerSetType()
164: @*/
165: PetscErrorCode VecTaggerSetUp(VecTagger tagger)
166: {

170:   if (tagger->setupcalled) return(0);
171:   if (!((PetscObject)tagger)->type_name) {VecTaggerSetType(tagger,VECTAGGERABSOLUTE);}
172:   if (tagger->ops->setup) {(*tagger->ops->setup)(tagger);}
173:   tagger->setupcalled = PETSC_TRUE;
174:   return(0);
175: }

177: /*@C
178:    VecTaggerSetFromOptions - set VecTagger options using the options database

180:    Logically Collective

182:    Input Arguments:
183: .  tagger - vec tagger

185:    Options Database Keys:
186: +  -vec_tagger_type       - implementation type, see VecTaggerSetType()
187: .  -vec_tagger_block_size - set the block size, see VecTaggerSetBlockSize()
188: -  -vec_tagger_invert     - invert the index set returned by VecTaggerComputeIS()

190:    Level: advanced

192: @*/
193: PetscErrorCode VecTaggerSetFromOptions(VecTagger tagger)
194: {
195:   VecTaggerType  deft;
196:   char           type[256];
198:   PetscBool      flg;

202:   PetscObjectOptionsBegin((PetscObject)tagger);
203:   deft = ((PetscObject)tagger)->type_name ? ((PetscObject)tagger)->type_name : VECTAGGERABSOLUTE;
204:   PetscOptionsFList("-vec_tagger_type","VecTagger implementation type","VecTaggerSetType",VecTaggerList,deft,type,256,&flg);
205:   VecTaggerSetType(tagger,flg ? type : deft);
206:   PetscOptionsInt("-vec_tagger_block_size","block size of the vectors the tagger operates on","VecTaggerSetBlockSize",tagger->blocksize,&tagger->blocksize,NULL);
207:   PetscOptionsBool("-vec_tagger_invert","invert the set of indices returned by VecTaggerComputeIS()","VecTaggerSetInvert",tagger->invert,&tagger->invert,NULL);
208:   if (tagger->ops->setfromoptions) {(*tagger->ops->setfromoptions)(PetscOptionsObject,tagger);}
209:   PetscOptionsEnd();
210:   return(0);
211: }

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

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

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

226:    Logically Collective

228:    Input Arguments:
229: +  tagger - vec tagger
230: -  blocksize - block size of the criteria used to tagger vectors

232:    Level: advanced

234: .seealso: VecTaggerComputeIS(), VecTaggerGetBlockSize(), VecSetBlockSize(), VecGetBlockSize()
235: @*/
236: PetscErrorCode VecTaggerSetBlockSize(VecTagger tagger, PetscInt blocksize)
237: {

242:   tagger->blocksize = blocksize;
243:   return(0);
244: }

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

249:    Logically Collective

251:    Input Arguments:
252: +  tagger - vec tagger
253: -  blocksize - block size of the vectors the tagger operates on

255:    Level: advanced

257: .seealso: VecTaggerComputeIS(), VecTaggerSetBlockSize(),
258: @*/
259: PetscErrorCode VecTaggerGetBlockSize(VecTagger tagger, PetscInt *blocksize)
260: {

265:   *blocksize = tagger->blocksize;
266:   return(0);
267: }

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

274:    Logically Collective

276:    Input Arguments:
277: +  tagger - vec tagger
278: -  invert - PETSC_TRUE to invert, PETSC_FALSE to use the indices as is

280:    Level: advanced

282: .seealso: VecTaggerComputeIS(), VecTaggerGetInvert()
283: @*/
284: PetscErrorCode VecTaggerSetInvert(VecTagger tagger, PetscBool invert)
285: {

290:   tagger->invert = invert;
291:   return(0);
292: }

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

297:    Logically Collective

299:    Input Arguments:
300: +  tagger - vec tagger
301: -  invert - PETSC_TRUE to invert, PETSC_FALSE to use the indices as is

303:    Level: advanced

305: .seealso: VecTaggerComputeIS(), VecTaggerSetInvert()
306: @*/
307: PetscErrorCode VecTaggerGetInvert(VecTagger tagger, PetscBool *invert)
308: {

313:   *invert = tagger->invert;
314:   return(0);
315: }

317: /*@C
318:    VecTaggerView - view a VecTagger context

320:    Collective

322:    Input Arguments:
323: +  tagger - vec tagger
324: -  viewer - viewer to display tagger, for example PETSC_VIEWER_STDOUT_WORLD

326:    Level: advanced

328: .seealso: VecTaggerCreate()
329: @*/
330: PetscErrorCode VecTaggerView(VecTagger tagger,PetscViewer viewer)
331: {
332:   PetscErrorCode    ierr;
333:   PetscBool         iascii;

337:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tagger),&viewer);}
340:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
341:   if (iascii) {
342:     PetscObjectPrintClassNamePrefixType((PetscObject)tagger,viewer);
343:     PetscViewerASCIIPushTab(viewer);
344:     PetscViewerASCIIPrintf(viewer,"Block size: %D\n",tagger->blocksize);
345:     if (tagger->ops->view) {(*tagger->ops->view)(tagger,viewer);}
346:     if (tagger->invert) {PetscViewerASCIIPrintf(viewer,"Inverting ISs.\n");}
347:     PetscViewerASCIIPopTab(viewer);
348:   }
349:   return(0);
350: }

352: /*@C
353:    VecTaggerComputeBoxes - If the tagged index set can be summarized as a list of boxes of values, returns that list

355:    Collective on VecTagger

357:    Input Aguments:
358: +  tagger - the VecTagger context
359: -  vec - the vec to tag

361:    Output Arguments:
362: +  numBoxes - the number of boxes in the tag definition
363: -  boxes - a newly allocated list of boxes.  This is a flat array of (BlockSize * numBoxes) pairs that the user can free with PetscFree().

365:    Notes:
366: .  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.

368:    Level: advanced

370: .seealso: VecTaggerComputeIS()
371: @*/
372: PetscErrorCode VecTaggerComputeBoxes(VecTagger tagger,Vec vec,PetscInt *numBoxes,VecTaggerBox **boxes)
373: {
374:   PetscInt       vls, tbs;

382:   VecGetLocalSize(vec,&vls);
383:   VecTaggerGetBlockSize(tagger,&tbs);
384:   if (vls % tbs) SETERRQ2(PetscObjectComm((PetscObject)tagger),PETSC_ERR_ARG_INCOMP,"vec local size %D is not a multiple of tagger block size %D",vls,tbs);
385:   if (tagger->ops->computeboxes) {(*tagger->ops->computeboxes) (tagger,vec,numBoxes,boxes);}
386:   else {
387:     const char *type;
388:     PetscObjectGetType ((PetscObject)tagger,&type);
389:     SETERRQ1(PetscObjectComm((PetscObject)tagger),PETSC_ERR_SUP,"VecTagger type %s does not compute value boxes",type);
390:   }
391:   return(0);
392: }

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

397:    Collective on VecTagger

399:    Input Aguments:
400: +  tagger - the VecTagger context
401: -  vec - the vec to tag

403:    Output Arguments:
404: .  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().

406:    Level: advanced

408: .seealso: VecTaggerComputeBoxes()
409: @*/
410: PetscErrorCode VecTaggerComputeIS(VecTagger tagger,Vec vec,IS *is)
411: {
412:   PetscInt       vls, tbs;

419:   VecGetLocalSize(vec,&vls);
420:   VecTaggerGetBlockSize(tagger,&tbs);
421:   if (vls % tbs) SETERRQ2(PetscObjectComm((PetscObject)tagger),PETSC_ERR_ARG_INCOMP,"vec local size %D is not a multiple of tagger block size %D",vls,tbs);
422:   if (tagger->ops->computeis) {(*tagger->ops->computeis) (tagger,vec,is);}
423:   else {
424:     SETERRQ(PetscObjectComm((PetscObject)tagger),PETSC_ERR_SUP,"VecTagger type does not compute ISs");
425:   }
426:   return(0);
427: }

429: PetscErrorCode VecTaggerComputeIS_FromBoxes(VecTagger tagger, Vec vec, IS *is)
430: { PetscInt       numBoxes;
431:   VecTaggerBox   *boxes;
432:   PetscInt       numTagged, offset;
433:   PetscInt       *tagged;
434:   PetscInt       bs, b, i, j, k, n;
435:   PetscBool      invert;
436:   const PetscScalar *vecArray;

440:   VecTaggerGetBlockSize(tagger,&bs);
441:   VecTaggerComputeBoxes(tagger,vec,&numBoxes,&boxes);
442:   VecGetArrayRead(vec, &vecArray);
443:   VecGetLocalSize(vec, &n);
444:   invert = tagger->invert;
445:   numTagged = 0;
446:   offset = 0;
447:   tagged = NULL;
448:   for (i = 0; i < 2; i++) {
449:     if (i) {
450:       PetscMalloc1(numTagged,&tagged);
451:     }
452:     for (j = 0; j < n; j++) {

454:       for (k = 0; k < numBoxes; k++) {
455:         for (b = 0; b < bs; b++) {
456:           PetscScalar  val = vecArray[j * bs + b];
457:           PetscInt     l = k * bs + b;
458:           VecTaggerBox box;
459:           PetscBool    in;

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