Actual source code: block.c

petsc-3.8.3 2017-12-09
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:     PetscViewerASCIIPushSynchronized(viewer);
159:     if (is->isperm) {
160:       PetscViewerASCIISynchronizedPrintf(viewer,"Block Index set is permutation\n");
161:     }
162:     PetscViewerASCIISynchronizedPrintf(viewer,"Block size %D\n",bs);
163:     PetscViewerASCIISynchronizedPrintf(viewer,"Number of block indices in set %D\n",n);
164:     PetscViewerASCIISynchronizedPrintf(viewer,"The first indices of each block are\n");
165:     for (i=0; i<n; i++) {
166:       PetscViewerASCIISynchronizedPrintf(viewer,"Block %D Index %D\n",i,idx[i]);
167:     }
168:     PetscViewerFlush(viewer);
169:     PetscViewerASCIIPopSynchronized(viewer);
170:   }
171:   return(0);
172: }

174: static PetscErrorCode ISSort_Block(IS is)
175: {
176:   IS_Block       *sub = (IS_Block*)is->data;
177:   PetscInt       bs, n;

181:   if (sub->sorted) return(0);
182:   PetscLayoutGetBlockSize(is->map, &bs);
183:   PetscLayoutGetLocalSize(is->map, &n);
184:   PetscSortInt(n/bs,sub->idx);
185:   sub->sorted = PETSC_TRUE;
186:   return(0);
187: }

189: static PetscErrorCode ISSortRemoveDups_Block(IS is)
190: {
191:   IS_Block       *sub = (IS_Block*)is->data;
192:   PetscInt       bs, n, nb;

196:   PetscLayoutGetBlockSize(is->map, &bs);
197:   PetscLayoutGetLocalSize(is->map, &n);
198:   nb   = n/bs;
199:   if (sub->sorted) {
200:     PetscSortedRemoveDupsInt(&nb,sub->idx);
201:   } else {
202:     PetscSortRemoveDupsInt(&nb,sub->idx);
203:   }
204:   PetscLayoutSetLocalSize(is->map, nb*bs);
205:   PetscLayoutSetSize(is->map, PETSC_DECIDE);
206:   PetscLayoutSetUp(is->map);
207:   sub->sorted = PETSC_TRUE;
208:   return(0);
209: }

211: static PetscErrorCode ISSorted_Block(IS is,PetscBool  *flg)
212: {
213:   IS_Block *sub = (IS_Block*)is->data;

216:   *flg = sub->sorted;
217:   return(0);
218: }

220: static PetscErrorCode ISDuplicate_Block(IS is,IS *newIS)
221: {
223:   IS_Block       *sub = (IS_Block*)is->data;
224:   PetscInt        bs, n;

227:   PetscLayoutGetBlockSize(is->map, &bs);
228:   PetscLayoutGetLocalSize(is->map, &n);
229:   n   /= bs;
230:   ISCreateBlock(PetscObjectComm((PetscObject)is),bs,n,sub->idx,PETSC_COPY_VALUES,newIS);
231:   return(0);
232: }

234: static PetscErrorCode ISIdentity_Block(IS is,PetscBool  *ident)
235: {
236:   IS_Block      *is_block = (IS_Block*)is->data;
237:   PetscInt       i,bs,n,*idx = is_block->idx;

241:   PetscLayoutGetBlockSize(is->map, &bs);
242:   PetscLayoutGetLocalSize(is->map, &n);
243:   n   /= bs;
244:   is->isidentity = PETSC_TRUE;
245:   *ident         = PETSC_TRUE;
246:   for (i=0; i<n; i++) {
247:     if (idx[i] != i) {
248:       is->isidentity = PETSC_FALSE;
249:       *ident         = PETSC_FALSE;
250:       return(0);
251:     }
252:   }
253:   return(0);
254: }

256: static PetscErrorCode ISCopy_Block(IS is,IS isy)
257: {
258:   IS_Block       *is_block = (IS_Block*)is->data,*isy_block = (IS_Block*)isy->data;
259:   PetscInt       bs, n, N, bsy, ny, Ny;

263:   PetscLayoutGetBlockSize(is->map, &bs);
264:   PetscLayoutGetLocalSize(is->map, &n);
265:   PetscLayoutGetSize(is->map, &N);
266:   PetscLayoutGetBlockSize(isy->map, &bsy);
267:   PetscLayoutGetLocalSize(isy->map, &ny);
268:   PetscLayoutGetSize(isy->map, &Ny);
269:   if (n != ny || N != Ny || bs != bsy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Index sets incompatible");
270:   isy_block->sorted = is_block->sorted;
271:   PetscMemcpy(isy_block->idx,is_block->idx,(n/bs)*sizeof(PetscInt));
272:   return(0);
273: }

275: static PetscErrorCode ISOnComm_Block(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
276: {
278:   IS_Block       *sub = (IS_Block*)is->data;
279:   PetscInt       bs, n;

282:   if (mode == PETSC_OWN_POINTER) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Cannot use PETSC_OWN_POINTER");
283:   PetscLayoutGetBlockSize(is->map, &bs);
284:   PetscLayoutGetLocalSize(is->map, &n);
285:   ISCreateBlock(comm,bs,n/bs,sub->idx,mode,newis);
286:   return(0);
287: }

289: static PetscErrorCode ISSetBlockSize_Block(IS is,PetscInt bs)
290: {

294:   PetscLayoutSetBlockSize(is->map, bs);
295:   return(0);
296: }

298: static PetscErrorCode ISToGeneral_Block(IS inis)
299: {
300:   IS_Block       *sub   = (IS_Block*)inis->data;
301:   PetscInt       bs,n;
302:   const PetscInt *idx;

306:   ISGetBlockSize(inis,&bs);
307:   ISGetLocalSize(inis,&n);
308:   ISGetIndices(inis,&idx);
309:   if (bs == 1) {
310:     PetscCopyMode mode = sub->allocated ? PETSC_OWN_POINTER : PETSC_USE_POINTER;
311:     sub->allocated = PETSC_FALSE; /* prevent deallocation when changing the subtype*/
312:     ISSetType(inis,ISGENERAL);
313:     ISGeneralSetIndices(inis,n,idx,mode);
314:   } else {
315:     ISSetType(inis,ISGENERAL);
316:     ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
317:   }
318:   return(0);
319: }


322: static struct _ISOps myops = { ISGetSize_Block,
323:                                ISGetLocalSize_Block,
324:                                ISGetIndices_Block,
325:                                ISRestoreIndices_Block,
326:                                ISInvertPermutation_Block,
327:                                ISSort_Block,
328:                                ISSortRemoveDups_Block,
329:                                ISSorted_Block,
330:                                ISDuplicate_Block,
331:                                ISDestroy_Block,
332:                                ISView_Block,
333:                                ISLoad_Default,
334:                                ISIdentity_Block,
335:                                ISCopy_Block,
336:                                ISToGeneral_Block,
337:                                ISOnComm_Block,
338:                                ISSetBlockSize_Block,
339:                                0,
340:                                ISLocate_Block};

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

345:    Collective on IS

347:    Input Parameters:
348: +  is - the index set
349: .  bs - number of elements in each block, one for each block and count of block not indices
350: .   n - the length of the index set (the number of blocks)
351: .  idx - the list of integers, these are by block, not by location
352: +  mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported


355:    Notes:
356:    When the communicator is not MPI_COMM_SELF, the operations on the
357:    index sets, IS, are NOT conceptually the same as MPI_Group operations.
358:    The index sets are then distributed sets of indices and thus certain operations
359:    on them are collective.

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

365:    Level: beginner

367:   Concepts: IS^block
368:   Concepts: index sets^block
369:   Concepts: block^index set

371: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
372: @*/
373: PetscErrorCode  ISBlockSetIndices(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
374: {

378:   PetscUseMethod(is,"ISBlockSetIndices_C",(IS,PetscInt,PetscInt,const PetscInt[],PetscCopyMode),(is,bs,n,idx,mode));
379:   return(0);
380: }

382: static PetscErrorCode  ISBlockSetIndices_Block(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
383: {
385:   PetscInt       i,min,max;
386:   IS_Block       *sub = (IS_Block*)is->data;

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

393:   PetscLayoutSetLocalSize(is->map, n*bs);
394:   PetscLayoutSetBlockSize(is->map, bs);
395:   PetscLayoutSetUp(is->map);

397:   if (sub->allocated) {PetscFree(sub->idx);}
398:   if (mode == PETSC_COPY_VALUES) {
399:     PetscMalloc1(n,&sub->idx);
400:     PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
401:     PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));
402:     sub->allocated = PETSC_TRUE;
403:   } else if (mode == PETSC_OWN_POINTER) {
404:     sub->idx = (PetscInt*) idx;
405:     PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
406:     sub->allocated = PETSC_TRUE;
407:   } else if (mode == PETSC_USE_POINTER) {
408:     sub->idx = (PetscInt*) idx;
409:     sub->allocated = PETSC_FALSE;
410:   }

412:   sub->sorted = PETSC_TRUE;
413:   for (i=1; i<n; i++) {
414:     if (idx[i] < idx[i-1]) {sub->sorted = PETSC_FALSE; break;}
415:   }
416:   if (n) {
417:     min = max = idx[0];
418:     for (i=1; i<n; i++) {
419:       if (idx[i] < min) min = idx[i];
420:       if (idx[i] > max) max = idx[i];
421:     }
422:     is->min = bs*min;
423:     is->max = bs*max+bs-1;
424:   } else {
425:     is->min = PETSC_MAX_INT;
426:     is->max = PETSC_MIN_INT;
427:   }
428:   is->isperm     = PETSC_FALSE;
429:   is->isidentity = PETSC_FALSE;
430:   return(0);
431: }

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

437:    Collective on MPI_Comm

439:    Input Parameters:
440: +  comm - the MPI communicator
441: .  bs - number of elements in each block
442: .  n - the length of the index set (the number of blocks)
443: .  idx - the list of integers, one for each block and count of block not indices
444: -  mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported in this routine

446:    Output Parameter:
447: .  is - the new index set

449:    Notes:
450:    When the communicator is not MPI_COMM_SELF, the operations on the
451:    index sets, IS, are NOT conceptually the same as MPI_Group operations.
452:    The index sets are then distributed sets of indices and thus certain operations
453:    on them are collective.

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

459:    Level: beginner

461:   Concepts: IS^block
462:   Concepts: index sets^block
463:   Concepts: block^index set

465: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
466: @*/
467: PetscErrorCode  ISCreateBlock(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
468: {

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

477:   ISCreate(comm,is);
478:   ISSetType(*is,ISBLOCK);
479:   ISBlockSetIndices(*is,bs,n,idx,mode);
480:   return(0);
481: }

483: static PetscErrorCode  ISBlockGetIndices_Block(IS is,const PetscInt *idx[])
484: {
485:   IS_Block *sub = (IS_Block*)is->data;

488:   *idx = sub->idx;
489:   return(0);
490: }

492: static PetscErrorCode  ISBlockRestoreIndices_Block(IS is,const PetscInt *idx[])
493: {
495:   return(0);
496: }

498: /*@C
499:    ISBlockGetIndices - Gets the indices associated with each block.

501:    Not Collective

503:    Input Parameter:
504: .  is - the index set

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

509:    Level: intermediate

511:    Concepts: IS^block
512:    Concepts: index sets^getting indices
513:    Concepts: index sets^block

515: .seealso: ISGetIndices(), ISBlockRestoreIndices()
516: @*/
517: PetscErrorCode  ISBlockGetIndices(IS is,const PetscInt *idx[])
518: {

522:   PetscUseMethod(is,"ISBlockGetIndices_C",(IS,const PetscInt*[]),(is,idx));
523:   return(0);
524: }

526: /*@C
527:    ISBlockRestoreIndices - Restores the indices associated with each block.

529:    Not Collective

531:    Input Parameter:
532: .  is - the index set

534:    Output Parameter:
535: .  idx - the integer indices

537:    Level: intermediate

539:    Concepts: IS^block
540:    Concepts: index sets^getting indices
541:    Concepts: index sets^block

543: .seealso: ISRestoreIndices(), ISBlockGetIndices()
544: @*/
545: PetscErrorCode  ISBlockRestoreIndices(IS is,const PetscInt *idx[])
546: {

550:   PetscUseMethod(is,"ISBlockRestoreIndices_C",(IS,const PetscInt*[]),(is,idx));
551:   return(0);
552: }

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

557:    Not Collective

559:    Input Parameter:
560: .  is - the index set

562:    Output Parameter:
563: .  size - the local number of blocks

565:    Level: intermediate

567:    Concepts: IS^block sizes
568:    Concepts: index sets^block sizes

570: .seealso: ISGetBlockSize(), ISBlockGetSize(), ISGetSize(), ISCreateBlock()
571: @*/
572: PetscErrorCode  ISBlockGetLocalSize(IS is,PetscInt *size)
573: {

577:   PetscUseMethod(is,"ISBlockGetLocalSize_C",(IS,PetscInt*),(is,size));
578:   return(0);
579: }

581: static PetscErrorCode  ISBlockGetLocalSize_Block(IS is,PetscInt *size)
582: {
583:   PetscInt       bs, n;

587:   PetscLayoutGetBlockSize(is->map, &bs);
588:   PetscLayoutGetLocalSize(is->map, &n);
589:   *size = n/bs;
590:   return(0);
591: }

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

596:    Not Collective

598:    Input Parameter:
599: .  is - the index set

601:    Output Parameter:
602: .  size - the global number of blocks

604:    Level: intermediate

606:    Concepts: IS^block sizes
607:    Concepts: index sets^block sizes

609: .seealso: ISGetBlockSize(), ISBlockGetLocalSize(), ISGetSize(), ISCreateBlock()
610: @*/
611: PetscErrorCode  ISBlockGetSize(IS is,PetscInt *size)
612: {

616:   PetscUseMethod(is,"ISBlockGetSize_C",(IS,PetscInt*),(is,size));
617:   return(0);
618: }

620: static PetscErrorCode  ISBlockGetSize_Block(IS is,PetscInt *size)
621: {
622:   PetscInt       bs, N;

626:   PetscLayoutGetBlockSize(is->map, &bs);
627:   PetscLayoutGetSize(is->map, &N);
628:   *size = N/bs;
629:   return(0);
630: }

632: PETSC_EXTERN PetscErrorCode ISCreate_Block(IS is)
633: {
635:   IS_Block       *sub;

638:   PetscNewLog(is,&sub);
639:   is->data = (void *) sub;
640:   PetscMemcpy(is->ops,&myops,sizeof(myops));
641:   PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",ISBlockSetIndices_Block);
642:   PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",ISBlockGetIndices_Block);
643:   PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",ISBlockRestoreIndices_Block);
644:   PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",ISBlockGetSize_Block);
645:   PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",ISBlockGetLocalSize_Block);
646:   return(0);
647: }