Actual source code: block.c

  1: /*$Id: block.c,v 1.54 2001/03/23 23:21:10 balay Exp $*/
  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 petscsys.h

  9: EXTERN int VecInitializePackage(char *);

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

 18: #undef __FUNCT__  
 20: int ISDestroy_Block(IS is)
 21: {
 22:   IS_Block *is_block = (IS_Block*)is->data;

 26:   PetscFree(is_block->idx);
 27:   PetscFree(is_block);
 28:   PetscLogObjectDestroy(is);
 29:   PetscHeaderDestroy(is); return(0);
 30: }

 32: #undef __FUNCT__  
 34: int ISGetIndices_Block(IS in,int **idx)
 35: {
 36:   IS_Block *sub = (IS_Block*)in->data;
 37:   int      i,j,k,bs = sub->bs,n = sub->n,*ii,*jj,ierr;

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

 56: #undef __FUNCT__  
 58: int ISRestoreIndices_Block(IS in,int **idx)
 59: {
 60:   IS_Block *sub = (IS_Block*)in->data;
 61:   int      ierr;

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

 74: #undef __FUNCT__  
 76: int ISGetSize_Block(IS is,int *size)
 77: {
 78:   IS_Block *sub = (IS_Block *)is->data;

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

 85: #undef __FUNCT__  
 87: int ISGetLocalSize_Block(IS is,int *size)
 88: {
 89:   IS_Block *sub = (IS_Block *)is->data;

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

 96: #undef __FUNCT__  
 98: int ISInvertPermutation_Block(IS is,int nlocal,IS *isout)
 99: {
100:   IS_Block *sub = (IS_Block *)is->data;
101:   int      i,ierr,*ii,n = sub->n,*idx = sub->idx,size;

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

119: #undef __FUNCT__  
121: int ISView_Block(IS is, PetscViewer viewer)
122: {
123:   IS_Block    *sub = (IS_Block *)is->data;
124:   int         i,n = sub->n,*idx = sub->idx,ierr;
125:   PetscTruth  isascii;

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

146: #undef __FUNCT__  
148: int ISSort_Block(IS is)
149: {
150:   IS_Block *sub = (IS_Block *)is->data;
151:   int      ierr;

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

160: #undef __FUNCT__  
162: int ISSorted_Block(IS is,PetscTruth *flg)
163: {
164:   IS_Block *sub = (IS_Block *)is->data;

167:   *flg = sub->sorted;
168:   return(0);
169: }

171: #undef __FUNCT__  
173: int ISDuplicate_Block(IS is,IS *newIS)
174: {
175:   int      ierr;
176:   IS_Block *sub = (IS_Block *)is->data;

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

183: #undef __FUNCT__  
185: int ISIdentity_Block(IS is,PetscTruth *ident)
186: {
187:   IS_Block *is_block = (IS_Block*)is->data;
188:   int      i,n = is_block->n,*idx = is_block->idx,bs = is_block->bs;

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

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

220:    Collective on MPI_Comm

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

228:    Output Parameter:
229: .  is - the new index set

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

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

241:    Level: beginner

243:   Concepts: IS^block
244:   Concepts: index sets^block
245:   Concepts: block^index set

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

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

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


292: #undef __FUNCT__  
294: /*@C
295:    ISBlockGetIndices - Gets the indices associated with each block.

297:    Not Collective

299:    Input Parameter:
300: .  is - the index set

302:    Output Parameter:
303: .  idx - the integer indices

305:    Level: intermediate

307:    Concepts: IS^block
308:    Concepts: index sets^getting indices
309:    Concepts: index sets^block

311: .seealso: ISGetIndices(), ISBlockRestoreIndices()
312: @*/
313: int ISBlockGetIndices(IS in,int *idx[])
314: {
315:   IS_Block *sub;

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

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

327: #undef __FUNCT__  
329: /*@C
330:    ISBlockRestoreIndices - Restores the indices associated with each block.

332:    Not Collective

334:    Input Parameter:
335: .  is - the index set

337:    Output Parameter:
338: .  idx - the integer indices

340:    Level: intermediate

342:    Concepts: IS^block
343:    Concepts: index sets^getting indices
344:    Concepts: index sets^block

346: .seealso: ISRestoreIndices(), ISBlockGetIndices()
347: @*/
348: int ISBlockRestoreIndices(IS is,int *idx[])
349: {
353:   if (is->type != IS_BLOCK) SETERRQ(PETSC_ERR_ARG_WRONG,"Not a block index set");
354:   return(0);
355: }

357: #undef __FUNCT__  
359: /*@
360:    ISBlockGetBlockSize - Returns the number of elements in a block.

362:    Not Collective

364:    Input Parameter:
365: .  is - the index set

367:    Output Parameter:
368: .  size - the number of elements in a block

370:    Level: intermediate

372:    Concepts: IS^block size
373:    Concepts: index sets^block size

375: .seealso: ISBlockGetSize(), ISGetSize(), ISBlock(), ISCreateBlock()
376: @*/
377: int ISBlockGetBlockSize(IS is,int *size)
378: {
379:   IS_Block *sub;

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

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

391: #undef __FUNCT__  
393: /*@C
394:    ISBlock - Checks whether an index set is blocked.

396:    Not Collective

398:    Input Parameter:
399: .  is - the index set

401:    Output Parameter:
402: .  flag - PETSC_TRUE if a block index set, else PETSC_FALSE

404:    Level: intermediate

406:    Concepts: IS^block
407:    Concepts: index sets^block

409: .seealso: ISBlockGetSize(), ISGetSize(), ISBlockGetBlockSize(), ISCreateBlock()
410: @*/
411: int ISBlock(IS is,PetscTruth *flag)
412: {
416:   if (is->type != IS_BLOCK) *flag = PETSC_FALSE;
417:   else                          *flag = PETSC_TRUE;
418:   return(0);
419: }

421: #undef __FUNCT__  
423: /*@
424:    ISBlockGetSize - Returns the number of blocks in the index set.

426:    Not Collective

428:    Input Parameter:
429: .  is - the index set

431:    Output Parameter:
432: .  size - the number of blocks

434:    Level: intermediate

436:    Concepts: IS^block sizes
437:    Concepts: index sets^block sizes

439: .seealso: ISBlockGetBlockSize(), ISGetSize(), ISBlock(), ISCreateBlock()
440: @*/
441: int ISBlockGetSize(IS is,int *size)
442: {
443:   IS_Block *sub;

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

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