Actual source code: block.c

petsc-master 2019-05-22
Report Typos and Errors

  2: /*
  3:      Provides the functions for index sets (IS) defined by a list of integers.
  4:    These are for blocks of data, each block is indicated with a single integer.
  5: */
  6:  #include <petsc/private/isimpl.h>
  7:  #include <petscviewer.h>

  9: typedef struct {
 10:   PetscBool sorted;      /* are the blocks sorted? */
 11:   PetscBool allocated;   /* did we allocate the index array ourselves? */
 12:   PetscInt  *idx;
 13: } IS_Block;

 15: static PetscErrorCode ISDestroy_Block(IS is)
 16: {
 17:   IS_Block       *sub = (IS_Block*)is->data;

 21:   if (sub->allocated) {PetscFree(sub->idx);}
 22:   PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",NULL);
 23:   PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",NULL);
 24:   PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",NULL);
 25:   PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",NULL);
 26:   PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",NULL);
 27:   PetscFree(is->data);
 28:   return(0);
 29: }

 31: static PetscErrorCode ISLocate_Block(IS is,PetscInt key,PetscInt *location)
 32: {
 33:   IS_Block       *sub = (IS_Block*)is->data;
 34:   PetscInt       numIdx, i, bs, bkey, mkey;

 38:   PetscLayoutGetBlockSize(is->map,&bs);
 39:   PetscLayoutGetSize(is->map,&numIdx);
 40:   numIdx /= bs;
 41:   bkey    = key / bs;
 42:   mkey    = key % bs;
 43:   if (mkey < 0) {
 44:     bkey--;
 45:     mkey += bs;
 46:   }
 47:   if (sub->sorted) {
 48:     PetscFindInt(bkey,numIdx,sub->idx,location);
 49:   } else {
 50:     const PetscInt *idx = sub->idx;

 52:     *location = -1;
 53:     for (i = 0; i < numIdx; i++) {
 54:       if (idx[i] == bkey) {
 55:         *location = i;
 56:         break;
 57:       }
 58:     }
 59:   }
 60:   if (*location >= 0) {
 61:     *location = *location * bs + mkey;
 62:   }
 63:   return(0);
 64: }

 66: static PetscErrorCode ISGetIndices_Block(IS in,const PetscInt *idx[])
 67: {
 68:   IS_Block       *sub = (IS_Block*)in->data;
 70:   PetscInt       i,j,k,bs,n,*ii,*jj;

 73:   PetscLayoutGetBlockSize(in->map, &bs);
 74:   PetscLayoutGetLocalSize(in->map, &n);
 75:   n   /= bs;
 76:   if (bs == 1) *idx = sub->idx;
 77:   else {
 78:     PetscMalloc1(bs*n,&jj);
 79:     *idx = jj;
 80:     k    = 0;
 81:     ii   = sub->idx;
 82:     for (i=0; i<n; i++)
 83:       for (j=0; j<bs; j++)
 84:         jj[k++] = bs*ii[i] + j;
 85:   }
 86:   return(0);
 87: }

 89: static PetscErrorCode ISRestoreIndices_Block(IS is,const PetscInt *idx[])
 90: {
 91:   IS_Block       *sub = (IS_Block*)is->data;
 92:   PetscInt       bs;

 96:   PetscLayoutGetBlockSize(is->map, &bs);
 97:   if (bs != 1) {
 98:     PetscFree(*(void**)idx);
 99:   } else {
100:     if (*idx != sub->idx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
101:   }
102:   return(0);
103: }

105: static PetscErrorCode ISGetSize_Block(IS is,PetscInt *size)
106: {

110:   PetscLayoutGetSize(is->map, size);
111:   return(0);
112: }

114: static PetscErrorCode ISGetLocalSize_Block(IS is,PetscInt *size)
115: {

119:   PetscLayoutGetLocalSize(is->map, size);
120:   return(0);
121: }

123: static PetscErrorCode ISInvertPermutation_Block(IS is,PetscInt nlocal,IS *isout)
124: {
125:   IS_Block       *sub = (IS_Block*)is->data;
126:   PetscInt       i,*ii,bs,n,*idx = sub->idx;
127:   PetscMPIInt    size;

131:   MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
132:   PetscLayoutGetBlockSize(is->map, &bs);
133:   PetscLayoutGetLocalSize(is->map, &n);
134:   n   /= bs;
135:   if (size == 1) {
136:     PetscMalloc1(n,&ii);
137:     for (i=0; i<n; i++) ii[idx[i]] = i;
138:     ISCreateBlock(PETSC_COMM_SELF,bs,n,ii,PETSC_OWN_POINTER,isout);
139:     ISSetPermutation(*isout);
140:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No inversion written yet for block IS");
141:   return(0);
142: }

144: static PetscErrorCode ISView_Block(IS is, PetscViewer viewer)
145: {
146:   IS_Block       *sub = (IS_Block*)is->data;
148:   PetscInt       i,bs,n,*idx = sub->idx;
149:   PetscBool      iascii;

152:   PetscLayoutGetBlockSize(is->map, &bs);
153:   PetscLayoutGetLocalSize(is->map, &n);
154:   n   /= bs;
155:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
156:   if (iascii) {
157:     PetscViewerFormat fmt;

159:     PetscViewerGetFormat(viewer,&fmt);
160:     if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
161:       IS             ist;
162:       const char     *name;
163:       const PetscInt *idx;
164:       PetscInt       n;

166:       PetscObjectGetName((PetscObject)is,&name);
167:       ISGetLocalSize(is,&n);
168:       ISGetIndices(is,&idx);
169:       ISCreateGeneral(PetscObjectComm((PetscObject)is),n,idx,PETSC_USE_POINTER,&ist);
170:       PetscObjectSetName((PetscObject)ist,name);
171:       ISView(ist,viewer);
172:       ISDestroy(&ist);
173:       ISRestoreIndices(is,&idx);
174:     } else {
175:       PetscViewerASCIIPushSynchronized(viewer);
176:       if (is->isperm) {
177:         PetscViewerASCIISynchronizedPrintf(viewer,"Block Index set is permutation\n");
178:       }
179:       PetscViewerASCIISynchronizedPrintf(viewer,"Block size %D\n",bs);
180:       PetscViewerASCIISynchronizedPrintf(viewer,"Number of block indices in set %D\n",n);
181:       PetscViewerASCIISynchronizedPrintf(viewer,"The first indices of each block are\n");
182:       for (i=0; i<n; i++) {
183:         PetscViewerASCIISynchronizedPrintf(viewer,"Block %D Index %D\n",i,idx[i]);
184:       }
185:       PetscViewerFlush(viewer);
186:       PetscViewerASCIIPopSynchronized(viewer);
187:     }
188:   }
189:   return(0);
190: }

192: static PetscErrorCode ISSort_Block(IS is)
193: {
194:   IS_Block       *sub = (IS_Block*)is->data;
195:   PetscInt       bs, n;

199:   if (sub->sorted) return(0);
200:   PetscLayoutGetBlockSize(is->map, &bs);
201:   PetscLayoutGetLocalSize(is->map, &n);
202:   PetscSortInt(n/bs,sub->idx);
203:   sub->sorted = PETSC_TRUE;
204:   return(0);
205: }

207: static PetscErrorCode ISSortRemoveDups_Block(IS is)
208: {
209:   IS_Block       *sub = (IS_Block*)is->data;
210:   PetscInt       bs, n, nb;

214:   PetscLayoutGetBlockSize(is->map, &bs);
215:   PetscLayoutGetLocalSize(is->map, &n);
216:   nb   = n/bs;
217:   if (sub->sorted) {
218:     PetscSortedRemoveDupsInt(&nb,sub->idx);
219:   } else {
220:     PetscSortRemoveDupsInt(&nb,sub->idx);
221:   }
222:   PetscLayoutSetLocalSize(is->map, nb*bs);
223:   PetscLayoutSetSize(is->map, PETSC_DECIDE);
224:   PetscLayoutSetUp(is->map);
225:   sub->sorted = PETSC_TRUE;
226:   return(0);
227: }

229: static PetscErrorCode ISSorted_Block(IS is,PetscBool  *flg)
230: {
231:   IS_Block *sub = (IS_Block*)is->data;

234:   *flg = sub->sorted;
235:   return(0);
236: }

238: static PetscErrorCode ISDuplicate_Block(IS is,IS *newIS)
239: {
241:   IS_Block       *sub = (IS_Block*)is->data;
242:   PetscInt        bs, n;

245:   PetscLayoutGetBlockSize(is->map, &bs);
246:   PetscLayoutGetLocalSize(is->map, &n);
247:   n   /= bs;
248:   ISCreateBlock(PetscObjectComm((PetscObject)is),bs,n,sub->idx,PETSC_COPY_VALUES,newIS);
249:   return(0);
250: }

252: static PetscErrorCode ISIdentity_Block(IS is,PetscBool  *ident)
253: {
254:   IS_Block      *is_block = (IS_Block*)is->data;
255:   PetscInt       i,bs,n,*idx = is_block->idx;

259:   PetscLayoutGetBlockSize(is->map, &bs);
260:   PetscLayoutGetLocalSize(is->map, &n);
261:   n   /= bs;
262:   is->isidentity = PETSC_TRUE;
263:   *ident         = PETSC_TRUE;
264:   for (i=0; i<n; i++) {
265:     if (idx[i] != i) {
266:       is->isidentity = PETSC_FALSE;
267:       *ident         = PETSC_FALSE;
268:       return(0);
269:     }
270:   }
271:   return(0);
272: }

274: static PetscErrorCode ISCopy_Block(IS is,IS isy)
275: {
276:   IS_Block       *is_block = (IS_Block*)is->data,*isy_block = (IS_Block*)isy->data;
277:   PetscInt       bs, n, N, bsy, ny, Ny;

281:   PetscLayoutGetBlockSize(is->map, &bs);
282:   PetscLayoutGetLocalSize(is->map, &n);
283:   PetscLayoutGetSize(is->map, &N);
284:   PetscLayoutGetBlockSize(isy->map, &bsy);
285:   PetscLayoutGetLocalSize(isy->map, &ny);
286:   PetscLayoutGetSize(isy->map, &Ny);
287:   if (n != ny || N != Ny || bs != bsy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Index sets incompatible");
288:   isy_block->sorted = is_block->sorted;
289:   PetscMemcpy(isy_block->idx,is_block->idx,(n/bs)*sizeof(PetscInt));
290:   return(0);
291: }

293: static PetscErrorCode ISOnComm_Block(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
294: {
296:   IS_Block       *sub = (IS_Block*)is->data;
297:   PetscInt       bs, n;

300:   if (mode == PETSC_OWN_POINTER) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Cannot use PETSC_OWN_POINTER");
301:   PetscLayoutGetBlockSize(is->map, &bs);
302:   PetscLayoutGetLocalSize(is->map, &n);
303:   ISCreateBlock(comm,bs,n/bs,sub->idx,mode,newis);
304:   return(0);
305: }

307: static PetscErrorCode ISSetBlockSize_Block(IS is,PetscInt bs)
308: {

312:   PetscLayoutSetBlockSize(is->map, bs);
313:   return(0);
314: }

316: static PetscErrorCode ISToGeneral_Block(IS inis)
317: {
318:   IS_Block       *sub   = (IS_Block*)inis->data;
319:   PetscInt       bs,n;
320:   const PetscInt *idx;

324:   ISGetBlockSize(inis,&bs);
325:   ISGetLocalSize(inis,&n);
326:   ISGetIndices(inis,&idx);
327:   if (bs == 1) {
328:     PetscCopyMode mode = sub->allocated ? PETSC_OWN_POINTER : PETSC_USE_POINTER;
329:     sub->allocated = PETSC_FALSE; /* prevent deallocation when changing the subtype*/
330:     ISSetType(inis,ISGENERAL);
331:     ISGeneralSetIndices(inis,n,idx,mode);
332:   } else {
333:     ISSetType(inis,ISGENERAL);
334:     ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
335:   }
336:   return(0);
337: }


340: static struct _ISOps myops = { ISGetSize_Block,
341:                                ISGetLocalSize_Block,
342:                                ISGetIndices_Block,
343:                                ISRestoreIndices_Block,
344:                                ISInvertPermutation_Block,
345:                                ISSort_Block,
346:                                ISSortRemoveDups_Block,
347:                                ISSorted_Block,
348:                                ISDuplicate_Block,
349:                                ISDestroy_Block,
350:                                ISView_Block,
351:                                ISLoad_Default,
352:                                ISIdentity_Block,
353:                                ISCopy_Block,
354:                                ISToGeneral_Block,
355:                                ISOnComm_Block,
356:                                ISSetBlockSize_Block,
357:                                0,
358:                                ISLocate_Block};

360: /*@
361:    ISBlockSetIndices - The indices are relative to entries, not blocks.

363:    Collective on IS

365:    Input Parameters:
366: +  is - the index set
367: .  bs - number of elements in each block, one for each block and count of block not indices
368: .   n - the length of the index set (the number of blocks)
369: .  idx - the list of integers, these are by block, not by location
370: +  mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported


373:    Notes:
374:    When the communicator is not MPI_COMM_SELF, the operations on the
375:    index sets, IS, are NOT conceptually the same as MPI_Group operations.
376:    The index sets are then distributed sets of indices and thus certain operations
377:    on them are collective.

379:    Example:
380:    If you wish to index the values {0,1,4,5}, then use
381:    a block size of 2 and idx of {0,2}.

383:    Level: beginner

385:   Concepts: IS^block
386:   Concepts: index sets^block
387:   Concepts: block^index set

389: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
390: @*/
391: PetscErrorCode  ISBlockSetIndices(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
392: {

396:   PetscUseMethod(is,"ISBlockSetIndices_C",(IS,PetscInt,PetscInt,const PetscInt[],PetscCopyMode),(is,bs,n,idx,mode));
397:   return(0);
398: }

400: static PetscErrorCode  ISBlockSetIndices_Block(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
401: {
403:   PetscInt       i,min,max;
404:   IS_Block       *sub = (IS_Block*)is->data;

407:   if (bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"block size < 1");
408:   if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");

411:   PetscLayoutSetLocalSize(is->map, n*bs);
412:   PetscLayoutSetBlockSize(is->map, bs);
413:   PetscLayoutSetUp(is->map);

415:   if (sub->allocated) {PetscFree(sub->idx);}
416:   if (mode == PETSC_COPY_VALUES) {
417:     PetscMalloc1(n,&sub->idx);
418:     PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
419:     PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));
420:     sub->allocated = PETSC_TRUE;
421:   } else if (mode == PETSC_OWN_POINTER) {
422:     sub->idx = (PetscInt*) idx;
423:     PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
424:     sub->allocated = PETSC_TRUE;
425:   } else if (mode == PETSC_USE_POINTER) {
426:     sub->idx = (PetscInt*) idx;
427:     sub->allocated = PETSC_FALSE;
428:   }

430:   sub->sorted = PETSC_TRUE;
431:   for (i=1; i<n; i++) {
432:     if (idx[i] < idx[i-1]) {sub->sorted = PETSC_FALSE; break;}
433:   }
434:   if (n) {
435:     min = max = idx[0];
436:     for (i=1; i<n; i++) {
437:       if (idx[i] < min) min = idx[i];
438:       if (idx[i] > max) max = idx[i];
439:     }
440:     is->min = bs*min;
441:     is->max = bs*max+bs-1;
442:   } else {
443:     is->min = PETSC_MAX_INT;
444:     is->max = PETSC_MIN_INT;
445:   }
446:   is->isperm     = PETSC_FALSE;
447:   is->isidentity = PETSC_FALSE;
448:   return(0);
449: }

451: /*@
452:    ISCreateBlock - Creates a data structure for an index set containing
453:    a list of integers. The indices are relative to entries, not blocks.

455:    Collective on MPI_Comm

457:    Input Parameters:
458: +  comm - the MPI communicator
459: .  bs - number of elements in each block
460: .  n - the length of the index set (the number of blocks)
461: .  idx - the list of integers, one for each block and count of block not indices
462: -  mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported in this routine

464:    Output Parameter:
465: .  is - the new index set

467:    Notes:
468:    When the communicator is not MPI_COMM_SELF, the operations on the
469:    index sets, IS, are NOT conceptually the same as MPI_Group operations.
470:    The index sets are then distributed sets of indices and thus certain operations
471:    on them are collective.

473:    Example:
474:    If you wish to index the values {0,1,6,7}, then use
475:    a block size of 2 and idx of {0,3}.

477:    Level: beginner

479:   Concepts: IS^block
480:   Concepts: index sets^block
481:   Concepts: block^index set

483: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
484: @*/
485: PetscErrorCode  ISCreateBlock(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
486: {

491:   if (bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"block size < 1");
492:   if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");

495:   ISCreate(comm,is);
496:   ISSetType(*is,ISBLOCK);
497:   ISBlockSetIndices(*is,bs,n,idx,mode);
498:   return(0);
499: }

501: static PetscErrorCode  ISBlockGetIndices_Block(IS is,const PetscInt *idx[])
502: {
503:   IS_Block *sub = (IS_Block*)is->data;

506:   *idx = sub->idx;
507:   return(0);
508: }

510: static PetscErrorCode  ISBlockRestoreIndices_Block(IS is,const PetscInt *idx[])
511: {
513:   return(0);
514: }

516: /*@C
517:    ISBlockGetIndices - Gets the indices associated with each block.

519:    Not Collective

521:    Input Parameter:
522: .  is - the index set

524:    Output Parameter:
525: .  idx - the integer indices, one for each block and count of block not indices

527:    Level: intermediate

529:    Concepts: IS^block
530:    Concepts: index sets^getting indices
531:    Concepts: index sets^block

533: .seealso: ISGetIndices(), ISBlockRestoreIndices()
534: @*/
535: PetscErrorCode  ISBlockGetIndices(IS is,const PetscInt *idx[])
536: {

540:   PetscUseMethod(is,"ISBlockGetIndices_C",(IS,const PetscInt*[]),(is,idx));
541:   return(0);
542: }

544: /*@C
545:    ISBlockRestoreIndices - Restores the indices associated with each block.

547:    Not Collective

549:    Input Parameter:
550: .  is - the index set

552:    Output Parameter:
553: .  idx - the integer indices

555:    Level: intermediate

557:    Concepts: IS^block
558:    Concepts: index sets^getting indices
559:    Concepts: index sets^block

561: .seealso: ISRestoreIndices(), ISBlockGetIndices()
562: @*/
563: PetscErrorCode  ISBlockRestoreIndices(IS is,const PetscInt *idx[])
564: {

568:   PetscUseMethod(is,"ISBlockRestoreIndices_C",(IS,const PetscInt*[]),(is,idx));
569:   return(0);
570: }

572: /*@
573:    ISBlockGetLocalSize - Returns the local number of blocks in the index set.

575:    Not Collective

577:    Input Parameter:
578: .  is - the index set

580:    Output Parameter:
581: .  size - the local number of blocks

583:    Level: intermediate

585:    Concepts: IS^block sizes
586:    Concepts: index sets^block sizes

588: .seealso: ISGetBlockSize(), ISBlockGetSize(), ISGetSize(), ISCreateBlock()
589: @*/
590: PetscErrorCode  ISBlockGetLocalSize(IS is,PetscInt *size)
591: {

595:   PetscUseMethod(is,"ISBlockGetLocalSize_C",(IS,PetscInt*),(is,size));
596:   return(0);
597: }

599: static PetscErrorCode  ISBlockGetLocalSize_Block(IS is,PetscInt *size)
600: {
601:   PetscInt       bs, n;

605:   PetscLayoutGetBlockSize(is->map, &bs);
606:   PetscLayoutGetLocalSize(is->map, &n);
607:   *size = n/bs;
608:   return(0);
609: }

611: /*@
612:    ISBlockGetSize - Returns the global number of blocks in the index set.

614:    Not Collective

616:    Input Parameter:
617: .  is - the index set

619:    Output Parameter:
620: .  size - the global number of blocks

622:    Level: intermediate

624:    Concepts: IS^block sizes
625:    Concepts: index sets^block sizes

627: .seealso: ISGetBlockSize(), ISBlockGetLocalSize(), ISGetSize(), ISCreateBlock()
628: @*/
629: PetscErrorCode  ISBlockGetSize(IS is,PetscInt *size)
630: {

634:   PetscUseMethod(is,"ISBlockGetSize_C",(IS,PetscInt*),(is,size));
635:   return(0);
636: }

638: static PetscErrorCode  ISBlockGetSize_Block(IS is,PetscInt *size)
639: {
640:   PetscInt       bs, N;

644:   PetscLayoutGetBlockSize(is->map, &bs);
645:   PetscLayoutGetSize(is->map, &N);
646:   *size = N/bs;
647:   return(0);
648: }

650: PETSC_EXTERN PetscErrorCode ISCreate_Block(IS is)
651: {
653:   IS_Block       *sub;

656:   PetscNewLog(is,&sub);
657:   is->data = (void *) sub;
658:   PetscMemcpy(is->ops,&myops,sizeof(myops));
659:   PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",ISBlockSetIndices_Block);
660:   PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",ISBlockGetIndices_Block);
661:   PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",ISBlockRestoreIndices_Block);
662:   PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",ISBlockGetSize_Block);
663:   PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",ISBlockGetLocalSize_Block);
664:   return(0);
665: }