Actual source code: block.c

  1: #define PETSCVEC_DLL
  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 src/vec/is/isimpl.h
 7:  #include petscvec.h

  9: typedef struct {
 10:   PetscInt        N,n;            /* number of blocks */
 11:   PetscTruth      sorted;       /* are the blocks sorted? */
 12:   PetscInt        *idx;
 13:   PetscInt        bs;           /* blocksize */
 14: } IS_Block;

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

 24:   PetscFree(is_block->idx);
 25:   PetscFree(is_block);
 26:   PetscHeaderDestroy(is);
 27:   return(0);
 28: }

 32: PetscErrorCode ISGetIndices_Block(IS in,PetscInt **idx)
 33: {
 34:   IS_Block       *sub = (IS_Block*)in->data;
 36:   PetscInt       i,j,k,bs = sub->bs,n = sub->n,*ii,*jj;

 39:   if (sub->bs == 1) {
 40:     *idx = sub->idx;
 41:   } else {
 42:     PetscMalloc(sub->bs*sub->n*sizeof(PetscInt),&jj);
 43:     *idx = jj;
 44:     k    = 0;
 45:     ii   = sub->idx;
 46:     for (i=0; i<n; i++) {
 47:       for (j=0; j<bs; j++) {
 48:         jj[k++] = ii[i] + j;
 49:       }
 50:     }
 51:   }
 52:   return(0);
 53: }

 57: PetscErrorCode ISRestoreIndices_Block(IS in,PetscInt **idx)
 58: {
 59:   IS_Block       *sub = (IS_Block*)in->data;

 63:   if (sub->bs != 1) {
 64:     PetscFree(*idx);
 65:   } else {
 66:     if (*idx !=  sub->idx) {
 67:       SETERRQ(PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
 68:     }
 69:   }
 70:   return(0);
 71: }

 75: PetscErrorCode ISGetSize_Block(IS is,PetscInt *size)
 76: {
 77:   IS_Block *sub = (IS_Block *)is->data;

 80:   *size = sub->bs*sub->N;
 81:   return(0);
 82: }

 86: PetscErrorCode ISGetLocalSize_Block(IS is,PetscInt *size)
 87: {
 88:   IS_Block *sub = (IS_Block *)is->data;

 91:   *size = sub->bs*sub->n;
 92:   return(0);
 93: }

 97: PetscErrorCode ISInvertPermutation_Block(IS is,PetscInt nlocal,IS *isout)
 98: {
 99:   IS_Block       *sub = (IS_Block *)is->data;
100:   PetscInt       i,*ii,n = sub->n,*idx = sub->idx;
101:   PetscMPIInt    size;

105:   MPI_Comm_size(is->comm,&size);
106:   if (size == 1) {
107:     PetscMalloc(n*sizeof(PetscInt),&ii);
108:     for (i=0; i<n; i++) {
109:       ii[idx[i]] = i;
110:     }
111:     ISCreateBlock(PETSC_COMM_SELF,sub->bs,n,ii,isout);
112:     ISSetPermutation(*isout);
113:     PetscFree(ii);
114:   } else {
115:     SETERRQ(PETSC_ERR_SUP,"No inversion written yet for block IS");
116:   }
117:   return(0);
118: }

122: PetscErrorCode ISView_Block(IS is, PetscViewer viewer)
123: {
124:   IS_Block       *sub = (IS_Block *)is->data;
126:   PetscInt       i,n = sub->n,*idx = sub->idx;
127:   PetscTruth     iascii;

130:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
131:   if (iascii) {
132:     if (is->isperm) {
133:       PetscViewerASCIISynchronizedPrintf(viewer,"Block Index set is permutation\n");
134:     }
135:     PetscViewerASCIISynchronizedPrintf(viewer,"Block size %D\n",sub->bs);
136:     PetscViewerASCIISynchronizedPrintf(viewer,"Number of block indices in set %D\n",n);
137:     PetscViewerASCIISynchronizedPrintf(viewer,"The first indices of each block are\n");
138:     for (i=0; i<n; i++) {
139:       PetscViewerASCIISynchronizedPrintf(viewer,"%D %D\n",i,idx[i]);
140:     }
141:     PetscViewerFlush(viewer);
142:   } else {
143:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for this object",((PetscObject)viewer)->type_name);
144:   }
145:   return(0);
146: }

150: PetscErrorCode ISSort_Block(IS is)
151: {
152:   IS_Block       *sub = (IS_Block *)is->data;

156:   if (sub->sorted) return(0);
157:   PetscSortInt(sub->n,sub->idx);
158:   sub->sorted = PETSC_TRUE;
159:   return(0);
160: }

164: PetscErrorCode ISSorted_Block(IS is,PetscTruth *flg)
165: {
166:   IS_Block *sub = (IS_Block *)is->data;

169:   *flg = sub->sorted;
170:   return(0);
171: }

175: PetscErrorCode ISDuplicate_Block(IS is,IS *newIS)
176: {
178:   IS_Block       *sub = (IS_Block *)is->data;

181:   ISCreateBlock(is->comm,sub->bs,sub->n,sub->idx,newIS);
182:   return(0);
183: }

187: PetscErrorCode ISIdentity_Block(IS is,PetscTruth *ident)
188: {
189:   IS_Block *is_block = (IS_Block*)is->data;
190:   PetscInt i,n = is_block->n,*idx = is_block->idx,bs = is_block->bs;

193:   is->isidentity = PETSC_TRUE;
194:   *ident         = PETSC_TRUE;
195:   for (i=0; i<n; i++) {
196:     if (idx[i] != bs*i) {
197:       is->isidentity = PETSC_FALSE;
198:       *ident         = PETSC_FALSE;
199:       return(0);
200:     }
201:   }
202:   return(0);
203: }

205: static struct _ISOps myops = { ISGetSize_Block,
206:                                ISGetLocalSize_Block,
207:                                ISGetIndices_Block,
208:                                ISRestoreIndices_Block,
209:                                ISInvertPermutation_Block,
210:                                ISSort_Block,
211:                                ISSorted_Block,
212:                                ISDuplicate_Block,
213:                                ISDestroy_Block,
214:                                ISView_Block,
215:                                ISIdentity_Block };
218: /*@
219:    ISCreateBlock - Creates a data structure for an index set containing
220:    a list of integers. The indices are relative to entries, not blocks. 

222:    Collective on MPI_Comm

224:    Input Parameters:
225: +  n - the length of the index set (the number of blocks)
226: .  bs - number of elements in each block
227: .  idx - the list of integers
228: -  comm - the MPI communicator

230:    Output Parameter:
231: .  is - the new index set

233:    Notes:
234:    When the communicator is not MPI_COMM_SELF, the operations on the 
235:    index sets, IS, are NOT conceptually the same as MPI_Group operations. 
236:    The index sets are then distributed sets of indices and thus certain operations
237:    on them are collective. 

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

243:    Level: beginner

245:   Concepts: IS^block
246:   Concepts: index sets^block
247:   Concepts: block^index set

249: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
250: @*/
251: PetscErrorCode  ISCreateBlock(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscInt idx[],IS *is)
252: {
254:   PetscInt       i,min,max;
255:   IS             Nindex;
256:   IS_Block       *sub;
257:   PetscTruth     sorted = PETSC_TRUE;

261:   if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
263:   *is = PETSC_NULL;
264: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
265:   VecInitializePackage(PETSC_NULL);
266: #endif

268:   PetscHeaderCreate(Nindex,_p_IS,struct _ISOps,IS_COOKIE,IS_BLOCK,"IS",comm,ISDestroy,ISView);
269:   PetscNew(IS_Block,&sub);
270:   PetscLogObjectMemory(Nindex,sizeof(IS_Block)+n*sizeof(PetscInt)+sizeof(struct _p_IS));
271:   PetscMalloc(n*sizeof(PetscInt),&sub->idx);
272:   sub->n = n;
273:   MPI_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,comm);
274:   for (i=1; i<n; i++) {
275:     if (idx[i] < idx[i-1]) {sorted = PETSC_FALSE; break;}
276:   }
277:   if (n) {min = max = idx[0];} else {min = max = 0;}
278:   for (i=1; i<n; i++) {
279:     if (idx[i] < min) min = idx[i];
280:     if (idx[i] > max) max = idx[i];
281:   }
282:   PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));
283:   sub->sorted     = sorted;
284:   sub->bs         = bs;
285:   Nindex->min     = min;
286:   Nindex->max     = max;
287:   Nindex->data    = (void*)sub;
288:   PetscMemcpy(Nindex->ops,&myops,sizeof(myops));
289:   Nindex->isperm  = PETSC_FALSE;
290:   *is = Nindex; return(0);
291: }


296: /*@C
297:    ISBlockGetIndices - Gets the indices associated with each block.

299:    Not Collective

301:    Input Parameter:
302: .  is - the index set

304:    Output Parameter:
305: .  idx - the integer indices

307:    Level: intermediate

309:    Concepts: IS^block
310:    Concepts: index sets^getting indices
311:    Concepts: index sets^block

313: .seealso: ISGetIndices(), ISBlockRestoreIndices()
314: @*/
315: PetscErrorCode  ISBlockGetIndices(IS in,PetscInt *idx[])
316: {
317:   IS_Block *sub;

322:   if (in->type != IS_BLOCK) SETERRQ(PETSC_ERR_ARG_WRONG,"Not a block index set");

324:   sub = (IS_Block*)in->data;
325:   *idx = sub->idx;
326:   return(0);
327: }

331: /*@C
332:    ISBlockRestoreIndices - Restores the indices associated with each block.

334:    Not Collective

336:    Input Parameter:
337: .  is - the index set

339:    Output Parameter:
340: .  idx - the integer indices

342:    Level: intermediate

344:    Concepts: IS^block
345:    Concepts: index sets^getting indices
346:    Concepts: index sets^block

348: .seealso: ISRestoreIndices(), ISBlockGetIndices()
349: @*/
350: PetscErrorCode  ISBlockRestoreIndices(IS is,PetscInt *idx[])
351: {
355:   if (is->type != IS_BLOCK) SETERRQ(PETSC_ERR_ARG_WRONG,"Not a block index set");
356:   return(0);
357: }

361: /*@
362:    ISBlockGetBlockSize - Returns the number of elements in a block.

364:    Not Collective

366:    Input Parameter:
367: .  is - the index set

369:    Output Parameter:
370: .  size - the number of elements in a block

372:    Level: intermediate

374:    Concepts: IS^block size
375:    Concepts: index sets^block size

377: .seealso: ISBlockGetSize(), ISGetSize(), ISBlock(), ISCreateBlock()
378: @*/
379: PetscErrorCode  ISBlockGetBlockSize(IS is,PetscInt *size)
380: {
381:   IS_Block *sub;

386:   if (is->type != IS_BLOCK) SETERRQ(PETSC_ERR_ARG_WRONG,"Not a block index set");

388:   sub = (IS_Block *)is->data;
389:   *size = sub->bs;
390:   return(0);
391: }

395: /*@
396:    ISBlock - Checks whether an index set is blocked.

398:    Not Collective

400:    Input Parameter:
401: .  is - the index set

403:    Output Parameter:
404: .  flag - PETSC_TRUE if a block index set, else PETSC_FALSE

406:    Level: intermediate

408:    Concepts: IS^block
409:    Concepts: index sets^block

411: .seealso: ISBlockGetSize(), ISGetSize(), ISBlockGetBlockSize(), ISCreateBlock()
412: @*/
413: PetscErrorCode  ISBlock(IS is,PetscTruth *flag)
414: {
418:   if (is->type != IS_BLOCK) *flag = PETSC_FALSE;
419:   else                          *flag = PETSC_TRUE;
420:   return(0);
421: }

425: /*@
426:    ISBlockGetSize - Returns the number of blocks in the index set.

428:    Not Collective

430:    Input Parameter:
431: .  is - the index set

433:    Output Parameter:
434: .  size - the number of blocks

436:    Level: intermediate

438:    Concepts: IS^block sizes
439:    Concepts: index sets^block sizes

441: .seealso: ISBlockGetBlockSize(), ISGetSize(), ISBlock(), ISCreateBlock()
442: @*/
443: PetscErrorCode  ISBlockGetSize(IS is,PetscInt *size)
444: {
445:   IS_Block *sub;

450:   if (is->type != IS_BLOCK) SETERRQ(PETSC_ERR_ARG_WRONG,"Not a block index set");

452:   sub = (IS_Block *)is->data;
453:   *size = sub->n;
454:   return(0);
455: }