Actual source code: block.c
petsc-3.8.4 2018-03-24
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>
7: #include <petscvec.h>
8: #include <petscviewer.h>
10: typedef struct {
11: PetscBool sorted; /* are the blocks sorted? */
12: PetscBool allocated; /* did we allocate the index array ourselves? */
13: PetscInt *idx;
14: } IS_Block;
16: static PetscErrorCode ISDestroy_Block(IS is)
17: {
18: IS_Block *sub = (IS_Block*)is->data;
22: if (sub->allocated) {PetscFree(sub->idx);}
23: PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",NULL);
24: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",NULL);
25: PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",NULL);
26: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",NULL);
27: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",NULL);
28: PetscFree(is->data);
29: return(0);
30: }
32: static PetscErrorCode ISLocate_Block(IS is,PetscInt key,PetscInt *location)
33: {
34: IS_Block *sub = (IS_Block*)is->data;
35: PetscInt numIdx, i, bs, bkey, mkey;
39: PetscLayoutGetBlockSize(is->map,&bs);
40: PetscLayoutGetSize(is->map,&numIdx);
41: numIdx /= bs;
42: bkey = key / bs;
43: mkey = key % bs;
44: if (mkey < 0) {
45: bkey--;
46: mkey += bs;
47: }
48: if (sub->sorted) {
49: PetscFindInt(bkey,numIdx,sub->idx,location);
50: } else {
51: const PetscInt *idx = sub->idx;
53: *location = -1;
54: for (i = 0; i < numIdx; i++) {
55: if (idx[i] == bkey) {
56: *location = i;
57: break;
58: }
59: }
60: }
61: if (*location >= 0) {
62: *location = *location * bs + mkey;
63: }
64: return(0);
65: }
67: static PetscErrorCode ISGetIndices_Block(IS in,const PetscInt *idx[])
68: {
69: IS_Block *sub = (IS_Block*)in->data;
71: PetscInt i,j,k,bs,n,*ii,*jj;
74: PetscLayoutGetBlockSize(in->map, &bs);
75: PetscLayoutGetLocalSize(in->map, &n);
76: n /= bs;
77: if (bs == 1) *idx = sub->idx;
78: else {
79: PetscMalloc1(bs*n,&jj);
80: *idx = jj;
81: k = 0;
82: ii = sub->idx;
83: for (i=0; i<n; i++)
84: for (j=0; j<bs; j++)
85: jj[k++] = bs*ii[i] + j;
86: }
87: return(0);
88: }
90: static PetscErrorCode ISRestoreIndices_Block(IS is,const PetscInt *idx[])
91: {
92: IS_Block *sub = (IS_Block*)is->data;
93: PetscInt bs;
97: PetscLayoutGetBlockSize(is->map, &bs);
98: if (bs != 1) {
99: PetscFree(*(void**)idx);
100: } else {
101: if (*idx != sub->idx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
102: }
103: return(0);
104: }
106: static PetscErrorCode ISGetSize_Block(IS is,PetscInt *size)
107: {
111: PetscLayoutGetSize(is->map, size);
112: return(0);
113: }
115: static PetscErrorCode ISGetLocalSize_Block(IS is,PetscInt *size)
116: {
120: PetscLayoutGetLocalSize(is->map, size);
121: return(0);
122: }
124: static PetscErrorCode ISInvertPermutation_Block(IS is,PetscInt nlocal,IS *isout)
125: {
126: IS_Block *sub = (IS_Block*)is->data;
127: PetscInt i,*ii,bs,n,*idx = sub->idx;
128: PetscMPIInt size;
132: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
133: PetscLayoutGetBlockSize(is->map, &bs);
134: PetscLayoutGetLocalSize(is->map, &n);
135: n /= bs;
136: if (size == 1) {
137: PetscMalloc1(n,&ii);
138: for (i=0; i<n; i++) ii[idx[i]] = i;
139: ISCreateBlock(PETSC_COMM_SELF,bs,n,ii,PETSC_OWN_POINTER,isout);
140: ISSetPermutation(*isout);
141: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No inversion written yet for block IS");
142: return(0);
143: }
145: static PetscErrorCode ISView_Block(IS is, PetscViewer viewer)
146: {
147: IS_Block *sub = (IS_Block*)is->data;
149: PetscInt i,bs,n,*idx = sub->idx;
150: PetscBool iascii;
153: PetscLayoutGetBlockSize(is->map, &bs);
154: PetscLayoutGetLocalSize(is->map, &n);
155: n /= bs;
156: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
157: if (iascii) {
158: PetscViewerASCIIPushSynchronized(viewer);
159: if (is->isperm) {
160: PetscViewerASCIISynchronizedPrintf(viewer,"Block Index set is permutation\n");
161: }
162: PetscViewerASCIISynchronizedPrintf(viewer,"Block size %D\n",bs);
163: PetscViewerASCIISynchronizedPrintf(viewer,"Number of block indices in set %D\n",n);
164: PetscViewerASCIISynchronizedPrintf(viewer,"The first indices of each block are\n");
165: for (i=0; i<n; i++) {
166: PetscViewerASCIISynchronizedPrintf(viewer,"Block %D Index %D\n",i,idx[i]);
167: }
168: PetscViewerFlush(viewer);
169: PetscViewerASCIIPopSynchronized(viewer);
170: }
171: return(0);
172: }
174: static PetscErrorCode ISSort_Block(IS is)
175: {
176: IS_Block *sub = (IS_Block*)is->data;
177: PetscInt bs, n;
181: if (sub->sorted) return(0);
182: PetscLayoutGetBlockSize(is->map, &bs);
183: PetscLayoutGetLocalSize(is->map, &n);
184: PetscSortInt(n/bs,sub->idx);
185: sub->sorted = PETSC_TRUE;
186: return(0);
187: }
189: static PetscErrorCode ISSortRemoveDups_Block(IS is)
190: {
191: IS_Block *sub = (IS_Block*)is->data;
192: PetscInt bs, n, nb;
196: PetscLayoutGetBlockSize(is->map, &bs);
197: PetscLayoutGetLocalSize(is->map, &n);
198: nb = n/bs;
199: if (sub->sorted) {
200: PetscSortedRemoveDupsInt(&nb,sub->idx);
201: } else {
202: PetscSortRemoveDupsInt(&nb,sub->idx);
203: }
204: PetscLayoutSetLocalSize(is->map, nb*bs);
205: PetscLayoutSetSize(is->map, PETSC_DECIDE);
206: PetscLayoutSetUp(is->map);
207: sub->sorted = PETSC_TRUE;
208: return(0);
209: }
211: static PetscErrorCode ISSorted_Block(IS is,PetscBool *flg)
212: {
213: IS_Block *sub = (IS_Block*)is->data;
216: *flg = sub->sorted;
217: return(0);
218: }
220: static PetscErrorCode ISDuplicate_Block(IS is,IS *newIS)
221: {
223: IS_Block *sub = (IS_Block*)is->data;
224: PetscInt bs, n;
227: PetscLayoutGetBlockSize(is->map, &bs);
228: PetscLayoutGetLocalSize(is->map, &n);
229: n /= bs;
230: ISCreateBlock(PetscObjectComm((PetscObject)is),bs,n,sub->idx,PETSC_COPY_VALUES,newIS);
231: return(0);
232: }
234: static PetscErrorCode ISIdentity_Block(IS is,PetscBool *ident)
235: {
236: IS_Block *is_block = (IS_Block*)is->data;
237: PetscInt i,bs,n,*idx = is_block->idx;
241: PetscLayoutGetBlockSize(is->map, &bs);
242: PetscLayoutGetLocalSize(is->map, &n);
243: n /= bs;
244: is->isidentity = PETSC_TRUE;
245: *ident = PETSC_TRUE;
246: for (i=0; i<n; i++) {
247: if (idx[i] != i) {
248: is->isidentity = PETSC_FALSE;
249: *ident = PETSC_FALSE;
250: return(0);
251: }
252: }
253: return(0);
254: }
256: static PetscErrorCode ISCopy_Block(IS is,IS isy)
257: {
258: IS_Block *is_block = (IS_Block*)is->data,*isy_block = (IS_Block*)isy->data;
259: PetscInt bs, n, N, bsy, ny, Ny;
263: PetscLayoutGetBlockSize(is->map, &bs);
264: PetscLayoutGetLocalSize(is->map, &n);
265: PetscLayoutGetSize(is->map, &N);
266: PetscLayoutGetBlockSize(isy->map, &bsy);
267: PetscLayoutGetLocalSize(isy->map, &ny);
268: PetscLayoutGetSize(isy->map, &Ny);
269: if (n != ny || N != Ny || bs != bsy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Index sets incompatible");
270: isy_block->sorted = is_block->sorted;
271: PetscMemcpy(isy_block->idx,is_block->idx,(n/bs)*sizeof(PetscInt));
272: return(0);
273: }
275: static PetscErrorCode ISOnComm_Block(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
276: {
278: IS_Block *sub = (IS_Block*)is->data;
279: PetscInt bs, n;
282: if (mode == PETSC_OWN_POINTER) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Cannot use PETSC_OWN_POINTER");
283: PetscLayoutGetBlockSize(is->map, &bs);
284: PetscLayoutGetLocalSize(is->map, &n);
285: ISCreateBlock(comm,bs,n/bs,sub->idx,mode,newis);
286: return(0);
287: }
289: static PetscErrorCode ISSetBlockSize_Block(IS is,PetscInt bs)
290: {
294: PetscLayoutSetBlockSize(is->map, bs);
295: return(0);
296: }
298: static PetscErrorCode ISToGeneral_Block(IS inis)
299: {
300: IS_Block *sub = (IS_Block*)inis->data;
301: PetscInt bs,n;
302: const PetscInt *idx;
306: ISGetBlockSize(inis,&bs);
307: ISGetLocalSize(inis,&n);
308: ISGetIndices(inis,&idx);
309: if (bs == 1) {
310: PetscCopyMode mode = sub->allocated ? PETSC_OWN_POINTER : PETSC_USE_POINTER;
311: sub->allocated = PETSC_FALSE; /* prevent deallocation when changing the subtype*/
312: ISSetType(inis,ISGENERAL);
313: ISGeneralSetIndices(inis,n,idx,mode);
314: } else {
315: ISSetType(inis,ISGENERAL);
316: ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
317: }
318: return(0);
319: }
322: static struct _ISOps myops = { ISGetSize_Block,
323: ISGetLocalSize_Block,
324: ISGetIndices_Block,
325: ISRestoreIndices_Block,
326: ISInvertPermutation_Block,
327: ISSort_Block,
328: ISSortRemoveDups_Block,
329: ISSorted_Block,
330: ISDuplicate_Block,
331: ISDestroy_Block,
332: ISView_Block,
333: ISLoad_Default,
334: ISIdentity_Block,
335: ISCopy_Block,
336: ISToGeneral_Block,
337: ISOnComm_Block,
338: ISSetBlockSize_Block,
339: 0,
340: ISLocate_Block};
342: /*@
343: ISBlockSetIndices - The indices are relative to entries, not blocks.
345: Collective on IS
347: Input Parameters:
348: + is - the index set
349: . bs - number of elements in each block, one for each block and count of block not indices
350: . n - the length of the index set (the number of blocks)
351: . idx - the list of integers, these are by block, not by location
352: + mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported
355: Notes:
356: When the communicator is not MPI_COMM_SELF, the operations on the
357: index sets, IS, are NOT conceptually the same as MPI_Group operations.
358: The index sets are then distributed sets of indices and thus certain operations
359: on them are collective.
361: Example:
362: If you wish to index the values {0,1,4,5}, then use
363: a block size of 2 and idx of {0,2}.
365: Level: beginner
367: Concepts: IS^block
368: Concepts: index sets^block
369: Concepts: block^index set
371: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
372: @*/
373: PetscErrorCode ISBlockSetIndices(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
374: {
378: PetscUseMethod(is,"ISBlockSetIndices_C",(IS,PetscInt,PetscInt,const PetscInt[],PetscCopyMode),(is,bs,n,idx,mode));
379: return(0);
380: }
382: static PetscErrorCode ISBlockSetIndices_Block(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
383: {
385: PetscInt i,min,max;
386: IS_Block *sub = (IS_Block*)is->data;
389: if (bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"block size < 1");
390: if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
393: PetscLayoutSetLocalSize(is->map, n*bs);
394: PetscLayoutSetBlockSize(is->map, bs);
395: PetscLayoutSetUp(is->map);
397: if (sub->allocated) {PetscFree(sub->idx);}
398: if (mode == PETSC_COPY_VALUES) {
399: PetscMalloc1(n,&sub->idx);
400: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
401: PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));
402: sub->allocated = PETSC_TRUE;
403: } else if (mode == PETSC_OWN_POINTER) {
404: sub->idx = (PetscInt*) idx;
405: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
406: sub->allocated = PETSC_TRUE;
407: } else if (mode == PETSC_USE_POINTER) {
408: sub->idx = (PetscInt*) idx;
409: sub->allocated = PETSC_FALSE;
410: }
412: sub->sorted = PETSC_TRUE;
413: for (i=1; i<n; i++) {
414: if (idx[i] < idx[i-1]) {sub->sorted = PETSC_FALSE; break;}
415: }
416: if (n) {
417: min = max = idx[0];
418: for (i=1; i<n; i++) {
419: if (idx[i] < min) min = idx[i];
420: if (idx[i] > max) max = idx[i];
421: }
422: is->min = bs*min;
423: is->max = bs*max+bs-1;
424: } else {
425: is->min = PETSC_MAX_INT;
426: is->max = PETSC_MIN_INT;
427: }
428: is->isperm = PETSC_FALSE;
429: is->isidentity = PETSC_FALSE;
430: return(0);
431: }
433: /*@
434: ISCreateBlock - Creates a data structure for an index set containing
435: a list of integers. The indices are relative to entries, not blocks.
437: Collective on MPI_Comm
439: Input Parameters:
440: + comm - the MPI communicator
441: . bs - number of elements in each block
442: . n - the length of the index set (the number of blocks)
443: . idx - the list of integers, one for each block and count of block not indices
444: - mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported in this routine
446: Output Parameter:
447: . is - the new index set
449: Notes:
450: When the communicator is not MPI_COMM_SELF, the operations on the
451: index sets, IS, are NOT conceptually the same as MPI_Group operations.
452: The index sets are then distributed sets of indices and thus certain operations
453: on them are collective.
455: Example:
456: If you wish to index the values {0,1,6,7}, then use
457: a block size of 2 and idx of {0,3}.
459: Level: beginner
461: Concepts: IS^block
462: Concepts: index sets^block
463: Concepts: block^index set
465: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
466: @*/
467: PetscErrorCode ISCreateBlock(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
468: {
473: if (bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"block size < 1");
474: if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
477: ISCreate(comm,is);
478: ISSetType(*is,ISBLOCK);
479: ISBlockSetIndices(*is,bs,n,idx,mode);
480: return(0);
481: }
483: static PetscErrorCode ISBlockGetIndices_Block(IS is,const PetscInt *idx[])
484: {
485: IS_Block *sub = (IS_Block*)is->data;
488: *idx = sub->idx;
489: return(0);
490: }
492: static PetscErrorCode ISBlockRestoreIndices_Block(IS is,const PetscInt *idx[])
493: {
495: return(0);
496: }
498: /*@C
499: ISBlockGetIndices - Gets the indices associated with each block.
501: Not Collective
503: Input Parameter:
504: . is - the index set
506: Output Parameter:
507: . idx - the integer indices, one for each block and count of block not indices
509: Level: intermediate
511: Concepts: IS^block
512: Concepts: index sets^getting indices
513: Concepts: index sets^block
515: .seealso: ISGetIndices(), ISBlockRestoreIndices()
516: @*/
517: PetscErrorCode ISBlockGetIndices(IS is,const PetscInt *idx[])
518: {
522: PetscUseMethod(is,"ISBlockGetIndices_C",(IS,const PetscInt*[]),(is,idx));
523: return(0);
524: }
526: /*@C
527: ISBlockRestoreIndices - Restores the indices associated with each block.
529: Not Collective
531: Input Parameter:
532: . is - the index set
534: Output Parameter:
535: . idx - the integer indices
537: Level: intermediate
539: Concepts: IS^block
540: Concepts: index sets^getting indices
541: Concepts: index sets^block
543: .seealso: ISRestoreIndices(), ISBlockGetIndices()
544: @*/
545: PetscErrorCode ISBlockRestoreIndices(IS is,const PetscInt *idx[])
546: {
550: PetscUseMethod(is,"ISBlockRestoreIndices_C",(IS,const PetscInt*[]),(is,idx));
551: return(0);
552: }
554: /*@
555: ISBlockGetLocalSize - Returns the local number of blocks in the index set.
557: Not Collective
559: Input Parameter:
560: . is - the index set
562: Output Parameter:
563: . size - the local number of blocks
565: Level: intermediate
567: Concepts: IS^block sizes
568: Concepts: index sets^block sizes
570: .seealso: ISGetBlockSize(), ISBlockGetSize(), ISGetSize(), ISCreateBlock()
571: @*/
572: PetscErrorCode ISBlockGetLocalSize(IS is,PetscInt *size)
573: {
577: PetscUseMethod(is,"ISBlockGetLocalSize_C",(IS,PetscInt*),(is,size));
578: return(0);
579: }
581: static PetscErrorCode ISBlockGetLocalSize_Block(IS is,PetscInt *size)
582: {
583: PetscInt bs, n;
587: PetscLayoutGetBlockSize(is->map, &bs);
588: PetscLayoutGetLocalSize(is->map, &n);
589: *size = n/bs;
590: return(0);
591: }
593: /*@
594: ISBlockGetSize - Returns the global number of blocks in the index set.
596: Not Collective
598: Input Parameter:
599: . is - the index set
601: Output Parameter:
602: . size - the global number of blocks
604: Level: intermediate
606: Concepts: IS^block sizes
607: Concepts: index sets^block sizes
609: .seealso: ISGetBlockSize(), ISBlockGetLocalSize(), ISGetSize(), ISCreateBlock()
610: @*/
611: PetscErrorCode ISBlockGetSize(IS is,PetscInt *size)
612: {
616: PetscUseMethod(is,"ISBlockGetSize_C",(IS,PetscInt*),(is,size));
617: return(0);
618: }
620: static PetscErrorCode ISBlockGetSize_Block(IS is,PetscInt *size)
621: {
622: PetscInt bs, N;
626: PetscLayoutGetBlockSize(is->map, &bs);
627: PetscLayoutGetSize(is->map, &N);
628: *size = N/bs;
629: return(0);
630: }
632: PETSC_EXTERN PetscErrorCode ISCreate_Block(IS is)
633: {
635: IS_Block *sub;
638: PetscNewLog(is,&sub);
639: is->data = (void *) sub;
640: PetscMemcpy(is->ops,&myops,sizeof(myops));
641: PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",ISBlockSetIndices_Block);
642: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",ISBlockGetIndices_Block);
643: PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",ISBlockRestoreIndices_Block);
644: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",ISBlockGetSize_Block);
645: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",ISBlockGetLocalSize_Block);
646: return(0);
647: }