Actual source code: block.c
petsc-3.6.4 2016-04-12
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: PetscBool sorted; /* are the blocks sorted? */
12: PetscBool borrowed_indices; /* do not free indices when IS is destroyed */
13: PetscInt *idx;
14: } IS_Block;
18: PetscErrorCode ISDestroy_Block(IS is)
19: {
20: IS_Block *is_block = (IS_Block*)is->data;
24: if (!is_block->borrowed_indices) {
25: PetscFree(is_block->idx);
26: }
27: PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",0);
28: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",0);
29: PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",0);
30: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",0);
31: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",0);
32: PetscFree(is->data);
33: return(0);
34: }
38: PetscErrorCode ISGetIndices_Block(IS in,const PetscInt *idx[])
39: {
40: IS_Block *sub = (IS_Block*)in->data;
42: PetscInt i,j,k,bs,n,*ii,*jj;
45: PetscLayoutGetBlockSize(in->map, &bs);
46: PetscLayoutGetLocalSize(in->map, &n);
47: n /= bs;
48: if (bs == 1) *idx = sub->idx;
49: else {
50: PetscMalloc1(bs*n,&jj);
51: *idx = jj;
52: k = 0;
53: ii = sub->idx;
54: for (i=0; i<n; i++)
55: for (j=0; j<bs; j++)
56: jj[k++] = bs*ii[i] + j;
57: }
58: return(0);
59: }
63: PetscErrorCode ISRestoreIndices_Block(IS is,const PetscInt *idx[])
64: {
65: IS_Block *sub = (IS_Block*)is->data;
66: PetscInt bs;
70: PetscLayoutGetBlockSize(is->map, &bs);
71: if (bs != 1) {
72: PetscFree(*(void**)idx);
73: } else {
74: if (*idx != sub->idx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
75: }
76: return(0);
77: }
81: PetscErrorCode ISGetSize_Block(IS is,PetscInt *size)
82: {
86: PetscLayoutGetSize(is->map, size);
87: return(0);
88: }
92: PetscErrorCode ISGetLocalSize_Block(IS is,PetscInt *size)
93: {
97: PetscLayoutGetLocalSize(is->map, size);
98: return(0);
99: }
103: PetscErrorCode ISInvertPermutation_Block(IS is,PetscInt nlocal,IS *isout)
104: {
105: IS_Block *sub = (IS_Block*)is->data;
106: PetscInt i,*ii,bs,n,*idx = sub->idx;
107: PetscMPIInt size;
111: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
112: PetscLayoutGetBlockSize(is->map, &bs);
113: PetscLayoutGetLocalSize(is->map, &n);
114: n /= bs;
115: if (size == 1) {
116: PetscMalloc1(n,&ii);
117: for (i=0; i<n; i++) ii[idx[i]] = i;
118: ISCreateBlock(PETSC_COMM_SELF,bs,n,ii,PETSC_OWN_POINTER,isout);
119: ISSetPermutation(*isout);
120: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No inversion written yet for block IS");
121: return(0);
122: }
126: PetscErrorCode ISView_Block(IS is, PetscViewer viewer)
127: {
128: IS_Block *sub = (IS_Block*)is->data;
130: PetscInt i,bs,n,*idx = sub->idx;
131: PetscBool iascii;
134: PetscLayoutGetBlockSize(is->map, &bs);
135: PetscLayoutGetLocalSize(is->map, &n);
136: n /= bs;
137: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
138: if (iascii) {
139: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
140: if (is->isperm) {
141: PetscViewerASCIISynchronizedPrintf(viewer,"Block Index set is permutation\n");
142: }
143: PetscViewerASCIISynchronizedPrintf(viewer,"Block size %D\n",bs);
144: PetscViewerASCIISynchronizedPrintf(viewer,"Number of block indices in set %D\n",n);
145: PetscViewerASCIISynchronizedPrintf(viewer,"The first indices of each block are\n");
146: for (i=0; i<n; i++) {
147: PetscViewerASCIISynchronizedPrintf(viewer,"Block %D Index %D\n",i,idx[i]);
148: }
149: PetscViewerFlush(viewer);
150: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
151: }
152: return(0);
153: }
157: PetscErrorCode ISSort_Block(IS is)
158: {
159: IS_Block *sub = (IS_Block*)is->data;
160: PetscInt bs, n;
164: if (sub->sorted) return(0);
165: PetscLayoutGetBlockSize(is->map, &bs);
166: PetscLayoutGetLocalSize(is->map, &n);
167: PetscSortInt(n/bs,sub->idx);
168: sub->sorted = PETSC_TRUE;
169: return(0);
170: }
174: PetscErrorCode ISSortRemoveDups_Block(IS is)
175: {
176: IS_Block *sub = (IS_Block*)is->data;
177: PetscInt bs, n, nb;
181: if (sub->sorted) return(0);
182: PetscLayoutGetBlockSize(is->map, &bs);
183: PetscLayoutGetLocalSize(is->map, &n);
184: nb = n/bs;
185: PetscSortRemoveDupsInt(&nb,sub->idx);
186: PetscLayoutSetLocalSize(is->map, nb*bs);
187: sub->sorted = PETSC_TRUE;
188: return(0);
189: }
193: PetscErrorCode ISSorted_Block(IS is,PetscBool *flg)
194: {
195: IS_Block *sub = (IS_Block*)is->data;
198: *flg = sub->sorted;
199: return(0);
200: }
204: PetscErrorCode ISDuplicate_Block(IS is,IS *newIS)
205: {
207: IS_Block *sub = (IS_Block*)is->data;
208: PetscInt bs, n;
211: PetscLayoutGetBlockSize(is->map, &bs);
212: PetscLayoutGetLocalSize(is->map, &n);
213: n /= bs;
214: ISCreateBlock(PetscObjectComm((PetscObject)is),bs,n,sub->idx,PETSC_COPY_VALUES,newIS);
215: return(0);
216: }
220: PetscErrorCode ISIdentity_Block(IS is,PetscBool *ident)
221: {
222: IS_Block *is_block = (IS_Block*)is->data;
223: PetscInt i,bs,n,*idx = is_block->idx;
227: PetscLayoutGetBlockSize(is->map, &bs);
228: PetscLayoutGetLocalSize(is->map, &n);
229: n /= bs;
230: is->isidentity = PETSC_TRUE;
231: *ident = PETSC_TRUE;
232: for (i=0; i<n; i++) {
233: if (idx[i] != i) {
234: is->isidentity = PETSC_FALSE;
235: *ident = PETSC_FALSE;
236: return(0);
237: }
238: }
239: return(0);
240: }
244: static PetscErrorCode ISCopy_Block(IS is,IS isy)
245: {
246: IS_Block *is_block = (IS_Block*)is->data,*isy_block = (IS_Block*)isy->data;
247: PetscInt bs, n, N, bsy, ny, Ny;
251: PetscLayoutGetBlockSize(is->map, &bs);
252: PetscLayoutGetLocalSize(is->map, &n);
253: PetscLayoutGetSize(is->map, &N);
254: PetscLayoutGetBlockSize(isy->map, &bsy);
255: PetscLayoutGetLocalSize(isy->map, &ny);
256: PetscLayoutGetSize(isy->map, &Ny);
257: if (n != ny || N != Ny || bs != bsy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Index sets incompatible");
258: isy_block->sorted = is_block->sorted;
259: PetscMemcpy(isy_block->idx,is_block->idx,(n/bs)*sizeof(PetscInt));
260: return(0);
261: }
265: static PetscErrorCode ISOnComm_Block(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
266: {
268: IS_Block *sub = (IS_Block*)is->data;
269: PetscInt bs, n;
272: if (mode == PETSC_OWN_POINTER) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Cannot use PETSC_OWN_POINTER");
273: PetscLayoutGetBlockSize(is->map, &bs);
274: PetscLayoutGetLocalSize(is->map, &n);
275: ISCreateBlock(comm,bs,n/bs,sub->idx,mode,newis);
276: return(0);
277: }
281: static PetscErrorCode ISSetBlockSize_Block(IS is,PetscInt bs)
282: {
286: PetscLayoutSetBlockSize(is->map, bs);
287: return(0);
288: }
292: static PetscErrorCode ISToGeneral_Block(IS inis)
293: {
294: IS_Block *sub = (IS_Block*)inis->data;
295: PetscInt bs,n;
296: const PetscInt *idx;
300: ISGetBlockSize(inis,&bs);
301: ISGetLocalSize(inis,&n);
302: ISGetIndices(inis,&idx);
303: if (bs == 1) {
304: PetscCopyMode mode = sub->borrowed_indices ? PETSC_USE_POINTER : PETSC_OWN_POINTER;
305: sub->borrowed_indices = PETSC_TRUE; /* prevent deallocation when changing the subtype*/
306: ISSetType(inis,ISGENERAL);
307: ISGeneralSetIndices(inis,n,idx,mode);
308: } else {
309: ISSetType(inis,ISGENERAL);
310: ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
311: }
312: return(0);
313: }
316: static struct _ISOps myops = { ISGetSize_Block,
317: ISGetLocalSize_Block,
318: ISGetIndices_Block,
319: ISRestoreIndices_Block,
320: ISInvertPermutation_Block,
321: ISSort_Block,
322: ISSortRemoveDups_Block,
323: ISSorted_Block,
324: ISDuplicate_Block,
325: ISDestroy_Block,
326: ISView_Block,
327: ISLoad_Default,
328: ISIdentity_Block,
329: ISCopy_Block,
330: ISToGeneral_Block,
331: ISOnComm_Block,
332: ISSetBlockSize_Block,
333: 0};
337: /*@
338: ISBlockSetIndices - The indices are relative to entries, not blocks.
340: Collective on IS
342: Input Parameters:
343: + is - the index set
344: . bs - number of elements in each block, one for each block and count of block not indices
345: . n - the length of the index set (the number of blocks)
346: . idx - the list of integers, these are by block, not by location
347: + mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported
350: Notes:
351: When the communicator is not MPI_COMM_SELF, the operations on the
352: index sets, IS, are NOT conceptually the same as MPI_Group operations.
353: The index sets are then distributed sets of indices and thus certain operations
354: on them are collective.
356: Example:
357: If you wish to index the values {0,1,4,5}, then use
358: a block size of 2 and idx of {0,2}.
360: Level: beginner
362: Concepts: IS^block
363: Concepts: index sets^block
364: Concepts: block^index set
366: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
367: @*/
368: PetscErrorCode ISBlockSetIndices(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
369: {
373: PetscUseMethod(is,"ISBlockSetIndices_C",(IS,PetscInt,PetscInt,const PetscInt[],PetscCopyMode),(is,bs,n,idx,mode));
374: return(0);
375: }
379: PetscErrorCode ISBlockSetIndices_Block(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
380: {
382: PetscInt i,min,max;
383: IS_Block *sub = (IS_Block*)is->data;
384: PetscBool sorted = PETSC_TRUE;
387: if (!sub->borrowed_indices) {
388: PetscFree(sub->idx);
389: } else {
390: sub->borrowed_indices = PETSC_FALSE;
391: }
392: PetscLayoutSetLocalSize(is->map, n*bs);
393: PetscLayoutSetBlockSize(is->map, bs);
394: PetscLayoutSetUp(is->map);
395: for (i=1; i<n; i++) {
396: if (idx[i] < idx[i-1]) {sorted = PETSC_FALSE; break;}
397: }
398: if (n) min = max = idx[0];
399: else min = max = 0;
400: for (i=1; i<n; i++) {
401: if (idx[i] < min) min = idx[i];
402: if (idx[i] > max) max = idx[i];
403: }
404: if (mode == PETSC_COPY_VALUES) {
405: PetscMalloc1(n,&sub->idx);
406: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
407: PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));
408: } else if (mode == PETSC_OWN_POINTER) {
409: sub->idx = (PetscInt*) idx;
410: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
411: } else if (mode == PETSC_USE_POINTER) {
412: sub->idx = (PetscInt*) idx;
413: sub->borrowed_indices = PETSC_TRUE;
414: }
416: sub->sorted = sorted;
417: is->min = bs*min;
418: is->max = bs*max+bs-1;
419: is->data = (void*)sub;
420: PetscMemcpy(is->ops,&myops,sizeof(myops));
421: is->isperm = PETSC_FALSE;
422: return(0);
423: }
427: /*@
428: ISCreateBlock - Creates a data structure for an index set containing
429: a list of integers. The indices are relative to entries, not blocks.
431: Collective on MPI_Comm
433: Input Parameters:
434: + comm - the MPI communicator
435: . bs - number of elements in each block
436: . n - the length of the index set (the number of blocks)
437: . idx - the list of integers, one for each block and count of block not indices
438: - mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported in this routine
440: Output Parameter:
441: . is - the new index set
443: Notes:
444: When the communicator is not MPI_COMM_SELF, the operations on the
445: index sets, IS, are NOT conceptually the same as MPI_Group operations.
446: The index sets are then distributed sets of indices and thus certain operations
447: on them are collective.
449: Example:
450: If you wish to index the values {0,1,6,7}, then use
451: a block size of 2 and idx of {0,3}.
453: Level: beginner
455: Concepts: IS^block
456: Concepts: index sets^block
457: Concepts: block^index set
459: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
460: @*/
461: PetscErrorCode ISCreateBlock(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
462: {
467: if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
470: ISCreate(comm,is);
471: ISSetType(*is,ISBLOCK);
472: ISBlockSetIndices(*is,bs,n,idx,mode);
473: return(0);
474: }
478: PetscErrorCode ISBlockGetIndices_Block(IS is,const PetscInt *idx[])
479: {
480: IS_Block *sub = (IS_Block*)is->data;
483: *idx = sub->idx;
484: return(0);
485: }
489: PetscErrorCode ISBlockRestoreIndices_Block(IS is,const PetscInt *idx[])
490: {
492: return(0);
493: }
497: /*@C
498: ISBlockGetIndices - Gets the indices associated with each block.
500: Not Collective
502: Input Parameter:
503: . is - the index set
505: Output Parameter:
506: . idx - the integer indices, one for each block and count of block not indices
508: Level: intermediate
510: Concepts: IS^block
511: Concepts: index sets^getting indices
512: Concepts: index sets^block
514: .seealso: ISGetIndices(), ISBlockRestoreIndices()
515: @*/
516: PetscErrorCode ISBlockGetIndices(IS is,const PetscInt *idx[])
517: {
521: PetscUseMethod(is,"ISBlockGetIndices_C",(IS,const PetscInt*[]),(is,idx));
522: return(0);
523: }
527: /*@C
528: ISBlockRestoreIndices - Restores the indices associated with each block.
530: Not Collective
532: Input Parameter:
533: . is - the index set
535: Output Parameter:
536: . idx - the integer indices
538: Level: intermediate
540: Concepts: IS^block
541: Concepts: index sets^getting indices
542: Concepts: index sets^block
544: .seealso: ISRestoreIndices(), ISBlockGetIndices()
545: @*/
546: PetscErrorCode ISBlockRestoreIndices(IS is,const PetscInt *idx[])
547: {
551: PetscUseMethod(is,"ISBlockRestoreIndices_C",(IS,const PetscInt*[]),(is,idx));
552: return(0);
553: }
557: /*@
558: ISBlockGetLocalSize - Returns the local number of blocks in the index set.
560: Not Collective
562: Input Parameter:
563: . is - the index set
565: Output Parameter:
566: . size - the local number of blocks
568: Level: intermediate
570: Concepts: IS^block sizes
571: Concepts: index sets^block sizes
573: .seealso: ISGetBlockSize(), ISBlockGetSize(), ISGetSize(), ISCreateBlock()
574: @*/
575: PetscErrorCode ISBlockGetLocalSize(IS is,PetscInt *size)
576: {
580: PetscUseMethod(is,"ISBlockGetLocalSize_C",(IS,PetscInt*),(is,size));
581: return(0);
582: }
586: PetscErrorCode ISBlockGetLocalSize_Block(IS is,PetscInt *size)
587: {
588: PetscInt bs, n;
592: PetscLayoutGetBlockSize(is->map, &bs);
593: PetscLayoutGetLocalSize(is->map, &n);
594: *size = n/bs;
595: return(0);
596: }
600: /*@
601: ISBlockGetSize - Returns the global number of blocks in the index set.
603: Not Collective
605: Input Parameter:
606: . is - the index set
608: Output Parameter:
609: . size - the global number of blocks
611: Level: intermediate
613: Concepts: IS^block sizes
614: Concepts: index sets^block sizes
616: .seealso: ISGetBlockSize(), ISBlockGetLocalSize(), ISGetSize(), ISCreateBlock()
617: @*/
618: PetscErrorCode ISBlockGetSize(IS is,PetscInt *size)
619: {
623: PetscUseMethod(is,"ISBlockGetSize_C",(IS,PetscInt*),(is,size));
624: return(0);
625: }
629: PetscErrorCode ISBlockGetSize_Block(IS is,PetscInt *size)
630: {
631: PetscInt bs, N;
635: PetscLayoutGetBlockSize(is->map, &bs);
636: PetscLayoutGetSize(is->map, &N);
637: *size = N/bs;
638: return(0);
639: }
643: PETSC_EXTERN PetscErrorCode ISCreate_Block(IS is)
644: {
646: IS_Block *sub;
649: PetscNewLog(is,&sub);
650: is->data = sub;
651: PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",ISBlockSetIndices_Block);
652: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",ISBlockGetIndices_Block);
653: PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",ISBlockRestoreIndices_Block);
654: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",ISBlockGetSize_Block);
655: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",ISBlockGetLocalSize_Block);
656: return(0);
657: }