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