Actual source code: block.c

petsc-3.5.1 2014-08-06
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>               /*I  "petscis.h"     I*/
  7: #include <petscvec.h>
  8: #include <petscviewer.h>

 10: typedef struct {
 11:   PetscBool sorted;             /* are the blocks sorted? */
 12:   PetscBool borrowed_indices;   /* do not free indices when IS is destroyed */
 13:   PetscInt  *idx;
 14: } IS_Block;

 18: PetscErrorCode ISDestroy_Block(IS is)
 19: {
 20:   IS_Block       *is_block = (IS_Block*)is->data;

 24:   if (!is_block->borrowed_indices) {
 25:     PetscFree(is_block->idx);
 26:   }
 27:   PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",0);
 28:   PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",0);
 29:   PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",0);
 30:   PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",0);
 31:   PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",0);
 32:   PetscFree(is->data);
 33:   return(0);
 34: }

 38: PetscErrorCode ISGetIndices_Block(IS in,const PetscInt *idx[])
 39: {
 40:   IS_Block       *sub = (IS_Block*)in->data;
 42:   PetscInt       i,j,k,bs,n,*ii,*jj;

 45:   PetscLayoutGetBlockSize(in->map, &bs);
 46:   PetscLayoutGetLocalSize(in->map, &n);
 47:   n   /= bs;
 48:   if (bs == 1) *idx = sub->idx;
 49:   else {
 50:     PetscMalloc1(bs*n,&jj);
 51:     *idx = jj;
 52:     k    = 0;
 53:     ii   = sub->idx;
 54:     for (i=0; i<n; i++)
 55:       for (j=0; j<bs; j++)
 56:         jj[k++] = bs*ii[i] + j;
 57:   }
 58:   return(0);
 59: }

 63: PetscErrorCode ISRestoreIndices_Block(IS is,const PetscInt *idx[])
 64: {
 65:   IS_Block       *sub = (IS_Block*)is->data;
 66:   PetscInt       bs;

 70:   PetscLayoutGetBlockSize(is->map, &bs);
 71:   if (bs != 1) {
 72:     PetscFree(*(void**)idx);
 73:   } else {
 74:     if (*idx != sub->idx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
 75:   }
 76:   return(0);
 77: }

 81: PetscErrorCode ISGetSize_Block(IS is,PetscInt *size)
 82: {

 86:   PetscLayoutGetSize(is->map, size);
 87:   return(0);
 88: }

 92: PetscErrorCode ISGetLocalSize_Block(IS is,PetscInt *size)
 93: {

 97:   PetscLayoutGetLocalSize(is->map, size);
 98:   return(0);
 99: }

103: PetscErrorCode ISInvertPermutation_Block(IS is,PetscInt nlocal,IS *isout)
104: {
105:   IS_Block       *sub = (IS_Block*)is->data;
106:   PetscInt       i,*ii,bs,n,*idx = sub->idx;
107:   PetscMPIInt    size;

111:   MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
112:   PetscLayoutGetBlockSize(is->map, &bs);
113:   PetscLayoutGetLocalSize(is->map, &n);
114:   n   /= bs;
115:   if (size == 1) {
116:     PetscMalloc1(n,&ii);
117:     for (i=0; i<n; i++) ii[idx[i]] = i;
118:     ISCreateBlock(PETSC_COMM_SELF,bs,n,ii,PETSC_OWN_POINTER,isout);
119:     ISSetPermutation(*isout);
120:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No inversion written yet for block IS");
121:   return(0);
122: }

126: PetscErrorCode ISView_Block(IS is, PetscViewer viewer)
127: {
128:   IS_Block       *sub = (IS_Block*)is->data;
130:   PetscInt       i,bs,n,*idx = sub->idx;
131:   PetscBool      iascii;

134:   PetscLayoutGetBlockSize(is->map, &bs);
135:   PetscLayoutGetLocalSize(is->map, &n);
136:   n   /= bs;
137:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
138:   if (iascii) {
139:     PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
140:     if (is->isperm) {
141:       PetscViewerASCIISynchronizedPrintf(viewer,"Block Index set is permutation\n");
142:     }
143:     PetscViewerASCIISynchronizedPrintf(viewer,"Block size %D\n",bs);
144:     PetscViewerASCIISynchronizedPrintf(viewer,"Number of block indices in set %D\n",n);
145:     PetscViewerASCIISynchronizedPrintf(viewer,"The first indices of each block are\n");
146:     for (i=0; i<n; i++) {
147:       PetscViewerASCIISynchronizedPrintf(viewer,"Block %D Index %D\n",i,idx[i]);
148:     }
149:     PetscViewerFlush(viewer);
150:     PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
151:   }
152:   return(0);
153: }

157: PetscErrorCode ISSort_Block(IS is)
158: {
159:   IS_Block       *sub = (IS_Block*)is->data;
160:   PetscInt       bs, n;

164:   if (sub->sorted) return(0);
165:   PetscLayoutGetBlockSize(is->map, &bs);
166:   PetscLayoutGetLocalSize(is->map, &n);
167:   PetscSortInt(n/bs,sub->idx);
168:   sub->sorted = PETSC_TRUE;
169:   return(0);
170: }

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

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

193: PetscErrorCode ISSorted_Block(IS is,PetscBool  *flg)
194: {
195:   IS_Block *sub = (IS_Block*)is->data;

198:   *flg = sub->sorted;
199:   return(0);
200: }

204: PetscErrorCode ISDuplicate_Block(IS is,IS *newIS)
205: {
207:   IS_Block       *sub = (IS_Block*)is->data;
208:   PetscInt        bs, n;

211:   PetscLayoutGetBlockSize(is->map, &bs);
212:   PetscLayoutGetLocalSize(is->map, &n);
213:   n   /= bs;
214:   ISCreateBlock(PetscObjectComm((PetscObject)is),bs,n,sub->idx,PETSC_COPY_VALUES,newIS);
215:   return(0);
216: }

220: PetscErrorCode ISIdentity_Block(IS is,PetscBool  *ident)
221: {
222:   IS_Block      *is_block = (IS_Block*)is->data;
223:   PetscInt       i,bs,n,*idx = is_block->idx;

227:   PetscLayoutGetBlockSize(is->map, &bs);
228:   PetscLayoutGetLocalSize(is->map, &n);
229:   n   /= bs;
230:   is->isidentity = PETSC_TRUE;
231:   *ident         = PETSC_TRUE;
232:   for (i=0; i<n; i++) {
233:     if (idx[i] != i) {
234:       is->isidentity = PETSC_FALSE;
235:       *ident         = PETSC_FALSE;
236:       return(0);
237:     }
238:   }
239:   return(0);
240: }

244: static PetscErrorCode ISCopy_Block(IS is,IS isy)
245: {
246:   IS_Block       *is_block = (IS_Block*)is->data,*isy_block = (IS_Block*)isy->data;
247:   PetscInt       bs, n, N, bsy, ny, Ny;

251:   PetscLayoutGetBlockSize(is->map, &bs);
252:   PetscLayoutGetLocalSize(is->map, &n);
253:   PetscLayoutGetSize(is->map, &N);
254:   PetscLayoutGetBlockSize(isy->map, &bsy);
255:   PetscLayoutGetLocalSize(isy->map, &ny);
256:   PetscLayoutGetSize(isy->map, &Ny);
257:   if (n != ny || N != Ny || bs != bsy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Index sets incompatible");
258:   isy_block->sorted = is_block->sorted;
259:   PetscMemcpy(isy_block->idx,is_block->idx,(n/bs)*sizeof(PetscInt));
260:   return(0);
261: }

265: static PetscErrorCode ISOnComm_Block(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
266: {
268:   IS_Block       *sub = (IS_Block*)is->data;
269:   PetscInt       bs, n;

272:   if (mode == PETSC_OWN_POINTER) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Cannot use PETSC_OWN_POINTER");
273:   PetscLayoutGetBlockSize(is->map, &bs);
274:   PetscLayoutGetLocalSize(is->map, &n);
275:   ISCreateBlock(comm,bs,n/bs,sub->idx,mode,newis);
276:   return(0);
277: }

281: static PetscErrorCode ISSetBlockSize_Block(IS is,PetscInt bs)
282: {

286:   PetscLayoutSetBlockSize(is->map, bs);
287:   return(0);
288: }

292: static PetscErrorCode ISToGeneral_Block(IS inis)
293: {
294:   IS_Block       *sub   = (IS_Block*)inis->data;
295:   PetscInt       bs,n;
296:   const PetscInt *idx;

300:   ISGetBlockSize(inis,&bs);
301:   ISGetLocalSize(inis,&n);
302:   ISGetIndices(inis,&idx);
303:   if (bs == 1) {
304:     PetscCopyMode mode = sub->borrowed_indices ? PETSC_USE_POINTER : PETSC_OWN_POINTER;
305:     sub->borrowed_indices = PETSC_TRUE; /* prevent deallocation when changing the subtype*/
306:     ISSetType(inis,ISGENERAL);
307:     ISGeneralSetIndices(inis,n,idx,mode);
308:   } else {
309:     ISSetType(inis,ISGENERAL);
310:     ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
311:   }
312:   return(0);
313: }


316: static struct _ISOps myops = { ISGetSize_Block,
317:                                ISGetLocalSize_Block,
318:                                ISGetIndices_Block,
319:                                ISRestoreIndices_Block,
320:                                ISInvertPermutation_Block,
321:                                ISSort_Block,
322:                                ISSortRemoveDups_Block,
323:                                ISSorted_Block,
324:                                ISDuplicate_Block,
325:                                ISDestroy_Block,
326:                                ISView_Block,
327:                                ISLoad_Default,
328:                                ISIdentity_Block,
329:                                ISCopy_Block,
330:                                ISToGeneral_Block,
331:                                ISOnComm_Block,
332:                                ISSetBlockSize_Block,
333:                                0};

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

340:    Collective on IS

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


350:    Notes:
351:    When the communicator is not MPI_COMM_SELF, the operations on the
352:    index sets, IS, are NOT conceptually the same as MPI_Group operations.
353:    The index sets are then distributed sets of indices and thus certain operations
354:    on them are collective.

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

360:    Level: beginner

362:   Concepts: IS^block
363:   Concepts: index sets^block
364:   Concepts: block^index set

366: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
367: @*/
368: PetscErrorCode  ISBlockSetIndices(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
369: {

373:   PetscUseMethod(is,"ISBlockSetIndices_C",(IS,PetscInt,PetscInt,const PetscInt[],PetscCopyMode),(is,bs,n,idx,mode));
374:   return(0);
375: }

379: PetscErrorCode  ISBlockSetIndices_Block(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
380: {
382:   PetscInt       i,min,max;
383:   IS_Block       *sub   = (IS_Block*)is->data;
384:   PetscBool      sorted = PETSC_TRUE;

387:   if (!sub->borrowed_indices) {
388:     PetscFree(sub->idx);
389:   } else {
390:     sub->borrowed_indices = PETSC_FALSE;
391:   }
392:   PetscLayoutSetLocalSize(is->map, n*bs);
393:   PetscLayoutSetBlockSize(is->map, bs);
394:   PetscLayoutSetUp(is->map);
395:   for (i=1; i<n; i++) {
396:     if (idx[i] < idx[i-1]) {sorted = PETSC_FALSE; break;}
397:   }
398:   if (n) min = max = idx[0];
399:   else   min = max = 0;
400:   for (i=1; i<n; i++) {
401:     if (idx[i] < min) min = idx[i];
402:     if (idx[i] > max) max = idx[i];
403:   }
404:   if (mode == PETSC_COPY_VALUES) {
405:     PetscMalloc1(n,&sub->idx);
406:     PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
407:     PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));
408:   } else if (mode == PETSC_OWN_POINTER) {
409:     sub->idx = (PetscInt*) idx;
410:     PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
411:   } else if (mode == PETSC_USE_POINTER) {
412:     sub->idx = (PetscInt*) idx;
413:     sub->borrowed_indices = PETSC_TRUE;
414:   }

416:   sub->sorted = sorted;
417:   is->min     = bs*min;
418:   is->max     = bs*max+bs-1;
419:   is->data    = (void*)sub;
420:   PetscMemcpy(is->ops,&myops,sizeof(myops));
421:   is->isperm  = PETSC_FALSE;
422:   return(0);
423: }

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

431:    Collective on MPI_Comm

433:    Input Parameters:
434: +  comm - the MPI communicator
435: .  bs - number of elements in each block
436: .  n - the length of the index set (the number of blocks)
437: .  idx - the list of integers, one for each block and count of block not indices
438: -  mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported in this routine

440:    Output Parameter:
441: .  is - the new index set

443:    Notes:
444:    When the communicator is not MPI_COMM_SELF, the operations on the
445:    index sets, IS, are NOT conceptually the same as MPI_Group operations.
446:    The index sets are then distributed sets of indices and thus certain operations
447:    on them are collective.

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

453:    Level: beginner

455:   Concepts: IS^block
456:   Concepts: index sets^block
457:   Concepts: block^index set

459: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
460: @*/
461: PetscErrorCode  ISCreateBlock(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
462: {

467:   if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");

470:   ISCreate(comm,is);
471:   ISSetType(*is,ISBLOCK);
472:   ISBlockSetIndices(*is,bs,n,idx,mode);
473:   return(0);
474: }

478: PetscErrorCode  ISBlockGetIndices_Block(IS is,const PetscInt *idx[])
479: {
480:   IS_Block *sub = (IS_Block*)is->data;

483:   *idx = sub->idx;
484:   return(0);
485: }

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

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

500:    Not Collective

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

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

508:    Level: intermediate

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

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

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

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

530:    Not Collective

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

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

538:    Level: intermediate

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

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

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

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

560:    Not Collective

562:    Input Parameter:
563: .  is - the index set

565:    Output Parameter:
566: .  size - the local number of blocks

568:    Level: intermediate

570:    Concepts: IS^block sizes
571:    Concepts: index sets^block sizes

573: .seealso: ISGetBlockSize(), ISBlockGetSize(), ISGetSize(), ISCreateBlock()
574: @*/
575: PetscErrorCode  ISBlockGetLocalSize(IS is,PetscInt *size)
576: {

580:   PetscUseMethod(is,"ISBlockGetLocalSize_C",(IS,PetscInt*),(is,size));
581:   return(0);
582: }

586: PetscErrorCode  ISBlockGetLocalSize_Block(IS is,PetscInt *size)
587: {
588:   PetscInt       bs, n;

592:   PetscLayoutGetBlockSize(is->map, &bs);
593:   PetscLayoutGetLocalSize(is->map, &n);
594:   *size = n/bs;
595:   return(0);
596: }

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

603:    Not Collective

605:    Input Parameter:
606: .  is - the index set

608:    Output Parameter:
609: .  size - the global number of blocks

611:    Level: intermediate

613:    Concepts: IS^block sizes
614:    Concepts: index sets^block sizes

616: .seealso: ISGetBlockSize(), ISBlockGetLocalSize(), ISGetSize(), ISCreateBlock()
617: @*/
618: PetscErrorCode  ISBlockGetSize(IS is,PetscInt *size)
619: {

623:   PetscUseMethod(is,"ISBlockGetSize_C",(IS,PetscInt*),(is,size));
624:   return(0);
625: }

629: PetscErrorCode  ISBlockGetSize_Block(IS is,PetscInt *size)
630: {
631:   PetscInt       bs, N;

635:   PetscLayoutGetBlockSize(is->map, &bs);
636:   PetscLayoutGetSize(is->map, &N);
637:   *size = N/bs;
638:   return(0);
639: }

643: PETSC_EXTERN PetscErrorCode ISCreate_Block(IS is)
644: {
646:   IS_Block       *sub;

649:   PetscNewLog(is,&sub);
650:   is->data = sub;
651:   PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",ISBlockSetIndices_Block);
652:   PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",ISBlockGetIndices_Block);
653:   PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",ISBlockRestoreIndices_Block);
654:   PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",ISBlockGetSize_Block);
655:   PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",ISBlockGetLocalSize_Block);
656:   return(0);
657: }