Actual source code: block.c

petsc-3.11.3 2019-06-26
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 <petscvec.h>
  8:  #include <petscviewer.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

364:    Collective on IS

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


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

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

384:    Level: beginner

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

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

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

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

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

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

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

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

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

456:    Collective on MPI_Comm

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

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

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

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

478:    Level: beginner

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

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

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

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

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

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

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

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

520:    Not Collective

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

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

528:    Level: intermediate

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

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

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

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

548:    Not Collective

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

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

556:    Level: intermediate

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

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

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

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

576:    Not Collective

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

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

584:    Level: intermediate

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

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

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

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

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

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

615:    Not Collective

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

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

623:    Level: intermediate

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

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

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

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

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

651: PETSC_EXTERN PetscErrorCode ISCreate_Block(IS is)
652: {
654:   IS_Block       *sub;

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