Actual source code: block.c
petsc-3.4.5 2014-06-29
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: PetscInt N,n; /* number of blocks */
12: PetscBool sorted; /* are the blocks sorted? */
13: PetscInt *idx;
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: PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",0);
26: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",0);
27: PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",0);
28: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",0);
29: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",0);
30: PetscFree(is->data);
31: return(0);
32: }
36: PetscErrorCode ISGetIndices_Block(IS in,const PetscInt *idx[])
37: {
38: IS_Block *sub = (IS_Block*)in->data;
40: PetscInt i,j,k,bs = in->bs,n = sub->n,*ii,*jj;
43: if (bs == 1) *idx = sub->idx;
44: else {
45: PetscMalloc(bs*n*sizeof(PetscInt),&jj);
46: *idx = jj;
47: k = 0;
48: ii = sub->idx;
49: for (i=0; i<n; i++)
50: for (j=0; j<bs; j++)
51: jj[k++] = bs*ii[i] + j;
52: }
53: return(0);
54: }
58: PetscErrorCode ISRestoreIndices_Block(IS in,const PetscInt *idx[])
59: {
60: IS_Block *sub = (IS_Block*)in->data;
64: if (in->bs != 1) {
65: PetscFree(*(void**)idx);
66: } else {
67: if (*idx != sub->idx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
68: }
69: return(0);
70: }
74: PetscErrorCode ISGetSize_Block(IS is,PetscInt *size)
75: {
76: IS_Block *sub = (IS_Block*)is->data;
79: *size = is->bs*sub->N;
80: return(0);
81: }
85: PetscErrorCode ISGetLocalSize_Block(IS is,PetscInt *size)
86: {
87: IS_Block *sub = (IS_Block*)is->data;
90: *size = is->bs*sub->n;
91: return(0);
92: }
96: PetscErrorCode ISInvertPermutation_Block(IS is,PetscInt nlocal,IS *isout)
97: {
98: IS_Block *sub = (IS_Block*)is->data;
99: PetscInt i,*ii,n = sub->n,*idx = sub->idx;
100: PetscMPIInt size;
104: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
105: if (size == 1) {
106: PetscMalloc(n*sizeof(PetscInt),&ii);
107: for (i=0; i<n; i++) ii[idx[i]] = i;
108: ISCreateBlock(PETSC_COMM_SELF,is->bs,n,ii,PETSC_OWN_POINTER,isout);
109: ISSetPermutation(*isout);
110: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No inversion written yet for block IS");
111: return(0);
112: }
116: PetscErrorCode ISView_Block(IS is, PetscViewer viewer)
117: {
118: IS_Block *sub = (IS_Block*)is->data;
120: PetscInt i,n = sub->n,*idx = sub->idx;
121: PetscBool iascii;
124: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
125: if (iascii) {
126: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
127: if (is->isperm) {
128: PetscViewerASCIISynchronizedPrintf(viewer,"Block Index set is permutation\n");
129: }
130: PetscViewerASCIISynchronizedPrintf(viewer,"Block size %D\n",is->bs);
131: PetscViewerASCIISynchronizedPrintf(viewer,"Number of block indices in set %D\n",n);
132: PetscViewerASCIISynchronizedPrintf(viewer,"The first indices of each block are\n");
133: for (i=0; i<n; i++) {
134: PetscViewerASCIISynchronizedPrintf(viewer,"Block %D Index %D\n",i,idx[i]);
135: }
136: PetscViewerFlush(viewer);
137: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
138: }
139: return(0);
140: }
144: PetscErrorCode ISSort_Block(IS is)
145: {
146: IS_Block *sub = (IS_Block*)is->data;
150: if (sub->sorted) return(0);
151: PetscSortInt(sub->n,sub->idx);
152: sub->sorted = PETSC_TRUE;
153: return(0);
154: }
158: PetscErrorCode ISSorted_Block(IS is,PetscBool *flg)
159: {
160: IS_Block *sub = (IS_Block*)is->data;
163: *flg = sub->sorted;
164: return(0);
165: }
169: PetscErrorCode ISDuplicate_Block(IS is,IS *newIS)
170: {
172: IS_Block *sub = (IS_Block*)is->data;
175: ISCreateBlock(PetscObjectComm((PetscObject)is),is->bs,sub->n,sub->idx,PETSC_COPY_VALUES,newIS);
176: return(0);
177: }
181: PetscErrorCode ISIdentity_Block(IS is,PetscBool *ident)
182: {
183: IS_Block *is_block = (IS_Block*)is->data;
184: PetscInt i,n = is_block->n,*idx = is_block->idx,bs = is->bs;
187: is->isidentity = PETSC_TRUE;
188: *ident = PETSC_TRUE;
189: for (i=0; i<n; i++) {
190: if (idx[i] != bs*i) {
191: is->isidentity = PETSC_FALSE;
192: *ident = PETSC_FALSE;
193: return(0);
194: }
195: }
196: return(0);
197: }
201: static PetscErrorCode ISCopy_Block(IS is,IS isy)
202: {
203: IS_Block *is_block = (IS_Block*)is->data,*isy_block = (IS_Block*)isy->data;
207: if (is_block->n != isy_block->n || is_block->N != isy_block->N || is->bs != isy->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Index sets incompatible");
208: isy_block->sorted = is_block->sorted;
209: PetscMemcpy(isy_block->idx,is_block->idx,is_block->n*sizeof(PetscInt));
210: return(0);
211: }
215: static PetscErrorCode ISOnComm_Block(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
216: {
218: IS_Block *sub = (IS_Block*)is->data;
221: if (mode == PETSC_OWN_POINTER) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Cannot use PETSC_OWN_POINTER");
222: ISCreateBlock(comm,is->bs,sub->n,sub->idx,mode,newis);
223: return(0);
224: }
228: static PetscErrorCode ISSetBlockSize_Block(IS is,PetscInt bs)
229: {
231: if (is->bs != bs) SETERRQ2(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_SIZ,"Cannot change block size for ISBLOCK from %D to %D",is->bs,bs);
232: return(0);
233: }
236: static struct _ISOps myops = { ISGetSize_Block,
237: ISGetLocalSize_Block,
238: ISGetIndices_Block,
239: ISRestoreIndices_Block,
240: ISInvertPermutation_Block,
241: ISSort_Block,
242: ISSorted_Block,
243: ISDuplicate_Block,
244: ISDestroy_Block,
245: ISView_Block,
246: ISIdentity_Block,
247: ISCopy_Block,
248: 0,
249: ISOnComm_Block,
250: ISSetBlockSize_Block,
251: 0};
255: /*@
256: ISBlockSetIndices - The indices are relative to entries, not blocks.
258: Collective on IS
260: Input Parameters:
261: + is - the index set
262: . bs - number of elements in each block, one for each block and count of block not indices
263: . n - the length of the index set (the number of blocks)
264: . idx - the list of integers, these are by block, not by location
265: + mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported
268: Notes:
269: When the communicator is not MPI_COMM_SELF, the operations on the
270: index sets, IS, are NOT conceptually the same as MPI_Group operations.
271: The index sets are then distributed sets of indices and thus certain operations
272: on them are collective.
274: Example:
275: If you wish to index the values {0,1,4,5}, then use
276: a block size of 2 and idx of {0,2}.
278: Level: beginner
280: Concepts: IS^block
281: Concepts: index sets^block
282: Concepts: block^index set
284: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
285: @*/
286: PetscErrorCode ISBlockSetIndices(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
287: {
291: PetscUseMethod(is,"ISBlockSetIndices_C",(IS,PetscInt,PetscInt,const PetscInt[],PetscCopyMode),(is,bs,n,idx,mode));
292: return(0);
293: }
297: PetscErrorCode ISBlockSetIndices_Block(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
298: {
300: PetscInt i,min,max;
301: IS_Block *sub = (IS_Block*)is->data;
302: PetscBool sorted = PETSC_TRUE;
305: PetscFree(sub->idx);
306: sub->n = n;
307: MPI_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)is));
308: for (i=1; i<n; i++) {
309: if (idx[i] < idx[i-1]) {sorted = PETSC_FALSE; break;}
310: }
311: if (n) min = max = idx[0];
312: else min = max = 0;
313: for (i=1; i<n; i++) {
314: if (idx[i] < min) min = idx[i];
315: if (idx[i] > max) max = idx[i];
316: }
317: if (mode == PETSC_COPY_VALUES) {
318: PetscMalloc(n*sizeof(PetscInt),&sub->idx);
319: PetscLogObjectMemory(is,n*sizeof(PetscInt));
320: PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));
321: } else if (mode == PETSC_OWN_POINTER) sub->idx = (PetscInt*) idx;
322: else SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Only supports PETSC_COPY_VALUES and PETSC_OWN_POINTER");
324: sub->sorted = sorted;
325: is->bs = bs;
326: is->min = bs*min;
327: is->max = bs*max+bs-1;
328: is->data = (void*)sub;
329: PetscMemcpy(is->ops,&myops,sizeof(myops));
330: is->isperm = PETSC_FALSE;
331: return(0);
332: }
336: /*@
337: ISCreateBlock - Creates a data structure for an index set containing
338: a list of integers. The indices are relative to entries, not blocks.
340: Collective on MPI_Comm
342: Input Parameters:
343: + comm - the MPI communicator
344: . bs - number of elements in each block
345: . n - the length of the index set (the number of blocks)
346: . idx - the list of integers, one for each block and count of block not indices
347: - mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported in this routine
349: Output Parameter:
350: . is - the new index set
352: Notes:
353: When the communicator is not MPI_COMM_SELF, the operations on the
354: index sets, IS, are NOT conceptually the same as MPI_Group operations.
355: The index sets are then distributed sets of indices and thus certain operations
356: on them are collective.
358: Example:
359: If you wish to index the values {0,1,6,7}, then use
360: a block size of 2 and idx of {0,3}.
362: Level: beginner
364: Concepts: IS^block
365: Concepts: index sets^block
366: Concepts: block^index set
368: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
369: @*/
370: PetscErrorCode ISCreateBlock(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
371: {
376: if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
379: ISCreate(comm,is);
380: ISSetType(*is,ISBLOCK);
381: ISBlockSetIndices(*is,bs,n,idx,mode);
382: return(0);
383: }
387: PetscErrorCode ISBlockGetIndices_Block(IS is,const PetscInt *idx[])
388: {
389: IS_Block *sub = (IS_Block*)is->data;
392: *idx = sub->idx;
393: return(0);
394: }
398: PetscErrorCode ISBlockRestoreIndices_Block(IS is,const PetscInt *idx[])
399: {
401: return(0);
402: }
406: /*@C
407: ISBlockGetIndices - Gets the indices associated with each block.
409: Not Collective
411: Input Parameter:
412: . is - the index set
414: Output Parameter:
415: . idx - the integer indices, one for each block and count of block not indices
417: Level: intermediate
419: Concepts: IS^block
420: Concepts: index sets^getting indices
421: Concepts: index sets^block
423: .seealso: ISGetIndices(), ISBlockRestoreIndices()
424: @*/
425: PetscErrorCode ISBlockGetIndices(IS is,const PetscInt *idx[])
426: {
430: PetscUseMethod(is,"ISBlockGetIndices_C",(IS,const PetscInt*[]),(is,idx));
431: return(0);
432: }
436: /*@C
437: ISBlockRestoreIndices - Restores the indices associated with each block.
439: Not Collective
441: Input Parameter:
442: . is - the index set
444: Output Parameter:
445: . idx - the integer indices
447: Level: intermediate
449: Concepts: IS^block
450: Concepts: index sets^getting indices
451: Concepts: index sets^block
453: .seealso: ISRestoreIndices(), ISBlockGetIndices()
454: @*/
455: PetscErrorCode ISBlockRestoreIndices(IS is,const PetscInt *idx[])
456: {
460: PetscUseMethod(is,"ISBlockRestoreIndices_C",(IS,const PetscInt*[]),(is,idx));
461: return(0);
462: }
466: /*@
467: ISBlockGetLocalSize - Returns the local number of blocks in the index set.
469: Not Collective
471: Input Parameter:
472: . is - the index set
474: Output Parameter:
475: . size - the local number of blocks
477: Level: intermediate
479: Concepts: IS^block sizes
480: Concepts: index sets^block sizes
482: .seealso: ISGetBlockSize(), ISBlockGetSize(), ISGetSize(), ISCreateBlock()
483: @*/
484: PetscErrorCode ISBlockGetLocalSize(IS is,PetscInt *size)
485: {
489: PetscUseMethod(is,"ISBlockGetLocalSize_C",(IS,PetscInt*),(is,size));
490: return(0);
491: }
495: PetscErrorCode ISBlockGetLocalSize_Block(IS is,PetscInt *size)
496: {
497: IS_Block *sub = (IS_Block*)is->data;
500: *size = sub->n;
501: return(0);
502: }
506: /*@
507: ISBlockGetSize - Returns the global number of blocks in the index set.
509: Not Collective
511: Input Parameter:
512: . is - the index set
514: Output Parameter:
515: . size - the global number of blocks
517: Level: intermediate
519: Concepts: IS^block sizes
520: Concepts: index sets^block sizes
522: .seealso: ISGetBlockSize(), ISBlockGetLocalSize(), ISGetSize(), ISCreateBlock()
523: @*/
524: PetscErrorCode ISBlockGetSize(IS is,PetscInt *size)
525: {
529: PetscUseMethod(is,"ISBlockGetSize_C",(IS,PetscInt*),(is,size));
530: return(0);
531: }
535: PetscErrorCode ISBlockGetSize_Block(IS is,PetscInt *size)
536: {
537: IS_Block *sub = (IS_Block*)is->data;
540: *size = sub->N;
541: return(0);
542: }
546: PetscErrorCode ISToGeneral_Block(IS inis)
547: {
549: const PetscInt *idx;
550: PetscInt n;
553: ISGetLocalSize(inis,&n);
554: ISGetIndices(inis,&idx);
555: ISSetType(inis,ISGENERAL);
556: ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
557: return(0);
558: }
562: PETSC_EXTERN PetscErrorCode ISCreate_Block(IS is)
563: {
565: IS_Block *sub;
568: PetscNewLog(is,IS_Block,&sub);
569: is->data = sub;
570: PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",ISBlockSetIndices_Block);
571: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",ISBlockGetIndices_Block);
572: PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",ISBlockRestoreIndices_Block);
573: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",ISBlockGetSize_Block);
574: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",ISBlockGetLocalSize_Block);
575: return(0);
576: }