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: }