Actual source code: block.c
petsc-3.10.5 2019-03-28
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: PetscViewerFormat fmt;
160: PetscViewerGetFormat(viewer,&fmt);
161: if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
162: IS ist;
163: const char *name;
164: const PetscInt *idx;
165: PetscInt n;
167: PetscObjectGetName((PetscObject)is,&name);
168: ISGetLocalSize(is,&n);
169: ISGetIndices(is,&idx);
170: ISCreateGeneral(PetscObjectComm((PetscObject)is),n,idx,PETSC_USE_POINTER,&ist);
171: PetscObjectSetName((PetscObject)ist,name);
172: ISView(ist,viewer);
173: ISDestroy(&ist);
174: ISRestoreIndices(is,&idx);
175: } else {
176: PetscViewerASCIIPushSynchronized(viewer);
177: if (is->isperm) {
178: PetscViewerASCIISynchronizedPrintf(viewer,"Block Index set is permutation\n");
179: }
180: PetscViewerASCIISynchronizedPrintf(viewer,"Block size %D\n",bs);
181: PetscViewerASCIISynchronizedPrintf(viewer,"Number of block indices in set %D\n",n);
182: PetscViewerASCIISynchronizedPrintf(viewer,"The first indices of each block are\n");
183: for (i=0; i<n; i++) {
184: PetscViewerASCIISynchronizedPrintf(viewer,"Block %D Index %D\n",i,idx[i]);
185: }
186: PetscViewerFlush(viewer);
187: PetscViewerASCIIPopSynchronized(viewer);
188: }
189: }
190: return(0);
191: }
193: static PetscErrorCode ISSort_Block(IS is)
194: {
195: IS_Block *sub = (IS_Block*)is->data;
196: PetscInt bs, n;
200: if (sub->sorted) return(0);
201: PetscLayoutGetBlockSize(is->map, &bs);
202: PetscLayoutGetLocalSize(is->map, &n);
203: PetscSortInt(n/bs,sub->idx);
204: sub->sorted = PETSC_TRUE;
205: return(0);
206: }
208: static PetscErrorCode ISSortRemoveDups_Block(IS is)
209: {
210: IS_Block *sub = (IS_Block*)is->data;
211: PetscInt bs, n, nb;
215: PetscLayoutGetBlockSize(is->map, &bs);
216: PetscLayoutGetLocalSize(is->map, &n);
217: nb = n/bs;
218: if (sub->sorted) {
219: PetscSortedRemoveDupsInt(&nb,sub->idx);
220: } else {
221: PetscSortRemoveDupsInt(&nb,sub->idx);
222: }
223: PetscLayoutSetLocalSize(is->map, nb*bs);
224: PetscLayoutSetSize(is->map, PETSC_DECIDE);
225: PetscLayoutSetUp(is->map);
226: sub->sorted = PETSC_TRUE;
227: return(0);
228: }
230: static PetscErrorCode ISSorted_Block(IS is,PetscBool *flg)
231: {
232: IS_Block *sub = (IS_Block*)is->data;
235: *flg = sub->sorted;
236: return(0);
237: }
239: static PetscErrorCode ISDuplicate_Block(IS is,IS *newIS)
240: {
242: IS_Block *sub = (IS_Block*)is->data;
243: PetscInt bs, n;
246: PetscLayoutGetBlockSize(is->map, &bs);
247: PetscLayoutGetLocalSize(is->map, &n);
248: n /= bs;
249: ISCreateBlock(PetscObjectComm((PetscObject)is),bs,n,sub->idx,PETSC_COPY_VALUES,newIS);
250: return(0);
251: }
253: static PetscErrorCode ISIdentity_Block(IS is,PetscBool *ident)
254: {
255: IS_Block *is_block = (IS_Block*)is->data;
256: PetscInt i,bs,n,*idx = is_block->idx;
260: PetscLayoutGetBlockSize(is->map, &bs);
261: PetscLayoutGetLocalSize(is->map, &n);
262: n /= bs;
263: is->isidentity = PETSC_TRUE;
264: *ident = PETSC_TRUE;
265: for (i=0; i<n; i++) {
266: if (idx[i] != i) {
267: is->isidentity = PETSC_FALSE;
268: *ident = PETSC_FALSE;
269: return(0);
270: }
271: }
272: return(0);
273: }
275: static PetscErrorCode ISCopy_Block(IS is,IS isy)
276: {
277: IS_Block *is_block = (IS_Block*)is->data,*isy_block = (IS_Block*)isy->data;
278: PetscInt bs, n, N, bsy, ny, Ny;
282: PetscLayoutGetBlockSize(is->map, &bs);
283: PetscLayoutGetLocalSize(is->map, &n);
284: PetscLayoutGetSize(is->map, &N);
285: PetscLayoutGetBlockSize(isy->map, &bsy);
286: PetscLayoutGetLocalSize(isy->map, &ny);
287: PetscLayoutGetSize(isy->map, &Ny);
288: if (n != ny || N != Ny || bs != bsy) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Index sets incompatible");
289: isy_block->sorted = is_block->sorted;
290: PetscMemcpy(isy_block->idx,is_block->idx,(n/bs)*sizeof(PetscInt));
291: return(0);
292: }
294: static PetscErrorCode ISOnComm_Block(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
295: {
297: IS_Block *sub = (IS_Block*)is->data;
298: PetscInt bs, n;
301: if (mode == PETSC_OWN_POINTER) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Cannot use PETSC_OWN_POINTER");
302: PetscLayoutGetBlockSize(is->map, &bs);
303: PetscLayoutGetLocalSize(is->map, &n);
304: ISCreateBlock(comm,bs,n/bs,sub->idx,mode,newis);
305: return(0);
306: }
308: static PetscErrorCode ISSetBlockSize_Block(IS is,PetscInt bs)
309: {
313: PetscLayoutSetBlockSize(is->map, bs);
314: return(0);
315: }
317: static PetscErrorCode ISToGeneral_Block(IS inis)
318: {
319: IS_Block *sub = (IS_Block*)inis->data;
320: PetscInt bs,n;
321: const PetscInt *idx;
325: ISGetBlockSize(inis,&bs);
326: ISGetLocalSize(inis,&n);
327: ISGetIndices(inis,&idx);
328: if (bs == 1) {
329: PetscCopyMode mode = sub->allocated ? PETSC_OWN_POINTER : PETSC_USE_POINTER;
330: sub->allocated = PETSC_FALSE; /* prevent deallocation when changing the subtype*/
331: ISSetType(inis,ISGENERAL);
332: ISGeneralSetIndices(inis,n,idx,mode);
333: } else {
334: ISSetType(inis,ISGENERAL);
335: ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
336: }
337: return(0);
338: }
341: static struct _ISOps myops = { ISGetSize_Block,
342: ISGetLocalSize_Block,
343: ISGetIndices_Block,
344: ISRestoreIndices_Block,
345: ISInvertPermutation_Block,
346: ISSort_Block,
347: ISSortRemoveDups_Block,
348: ISSorted_Block,
349: ISDuplicate_Block,
350: ISDestroy_Block,
351: ISView_Block,
352: ISLoad_Default,
353: ISIdentity_Block,
354: ISCopy_Block,
355: ISToGeneral_Block,
356: ISOnComm_Block,
357: ISSetBlockSize_Block,
358: 0,
359: ISLocate_Block};
361: /*@
362: ISBlockSetIndices - The indices are relative to entries, not blocks.
364: Collective on IS
366: Input Parameters:
367: + is - the index set
368: . bs - number of elements in each block, one for each block and count of block not indices
369: . n - the length of the index set (the number of blocks)
370: . idx - the list of integers, these are by block, not by location
371: + mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported
374: Notes:
375: When the communicator is not MPI_COMM_SELF, the operations on the
376: index sets, IS, are NOT conceptually the same as MPI_Group operations.
377: The index sets are then distributed sets of indices and thus certain operations
378: on them are collective.
380: Example:
381: If you wish to index the values {0,1,4,5}, then use
382: a block size of 2 and idx of {0,2}.
384: Level: beginner
386: Concepts: IS^block
387: Concepts: index sets^block
388: Concepts: block^index set
390: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
391: @*/
392: PetscErrorCode ISBlockSetIndices(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
393: {
397: PetscUseMethod(is,"ISBlockSetIndices_C",(IS,PetscInt,PetscInt,const PetscInt[],PetscCopyMode),(is,bs,n,idx,mode));
398: return(0);
399: }
401: static PetscErrorCode ISBlockSetIndices_Block(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
402: {
404: PetscInt i,min,max;
405: IS_Block *sub = (IS_Block*)is->data;
408: if (bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"block size < 1");
409: if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
412: PetscLayoutSetLocalSize(is->map, n*bs);
413: PetscLayoutSetBlockSize(is->map, bs);
414: PetscLayoutSetUp(is->map);
416: if (sub->allocated) {PetscFree(sub->idx);}
417: if (mode == PETSC_COPY_VALUES) {
418: PetscMalloc1(n,&sub->idx);
419: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
420: PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));
421: sub->allocated = PETSC_TRUE;
422: } else if (mode == PETSC_OWN_POINTER) {
423: sub->idx = (PetscInt*) idx;
424: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
425: sub->allocated = PETSC_TRUE;
426: } else if (mode == PETSC_USE_POINTER) {
427: sub->idx = (PetscInt*) idx;
428: sub->allocated = PETSC_FALSE;
429: }
431: sub->sorted = PETSC_TRUE;
432: for (i=1; i<n; i++) {
433: if (idx[i] < idx[i-1]) {sub->sorted = PETSC_FALSE; break;}
434: }
435: if (n) {
436: min = max = idx[0];
437: for (i=1; i<n; i++) {
438: if (idx[i] < min) min = idx[i];
439: if (idx[i] > max) max = idx[i];
440: }
441: is->min = bs*min;
442: is->max = bs*max+bs-1;
443: } else {
444: is->min = PETSC_MAX_INT;
445: is->max = PETSC_MIN_INT;
446: }
447: is->isperm = PETSC_FALSE;
448: is->isidentity = PETSC_FALSE;
449: return(0);
450: }
452: /*@
453: ISCreateBlock - Creates a data structure for an index set containing
454: a list of integers. The indices are relative to entries, not blocks.
456: Collective on MPI_Comm
458: Input Parameters:
459: + comm - the MPI communicator
460: . bs - number of elements in each block
461: . n - the length of the index set (the number of blocks)
462: . idx - the list of integers, one for each block and count of block not indices
463: - mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported in this routine
465: Output Parameter:
466: . is - the new index set
468: Notes:
469: When the communicator is not MPI_COMM_SELF, the operations on the
470: index sets, IS, are NOT conceptually the same as MPI_Group operations.
471: The index sets are then distributed sets of indices and thus certain operations
472: on them are collective.
474: Example:
475: If you wish to index the values {0,1,6,7}, then use
476: a block size of 2 and idx of {0,3}.
478: Level: beginner
480: Concepts: IS^block
481: Concepts: index sets^block
482: Concepts: block^index set
484: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
485: @*/
486: PetscErrorCode ISCreateBlock(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
487: {
492: if (bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"block size < 1");
493: if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
496: ISCreate(comm,is);
497: ISSetType(*is,ISBLOCK);
498: ISBlockSetIndices(*is,bs,n,idx,mode);
499: return(0);
500: }
502: static PetscErrorCode ISBlockGetIndices_Block(IS is,const PetscInt *idx[])
503: {
504: IS_Block *sub = (IS_Block*)is->data;
507: *idx = sub->idx;
508: return(0);
509: }
511: static PetscErrorCode ISBlockRestoreIndices_Block(IS is,const PetscInt *idx[])
512: {
514: return(0);
515: }
517: /*@C
518: ISBlockGetIndices - Gets the indices associated with each block.
520: Not Collective
522: Input Parameter:
523: . is - the index set
525: Output Parameter:
526: . idx - the integer indices, one for each block and count of block not indices
528: Level: intermediate
530: Concepts: IS^block
531: Concepts: index sets^getting indices
532: Concepts: index sets^block
534: .seealso: ISGetIndices(), ISBlockRestoreIndices()
535: @*/
536: PetscErrorCode ISBlockGetIndices(IS is,const PetscInt *idx[])
537: {
541: PetscUseMethod(is,"ISBlockGetIndices_C",(IS,const PetscInt*[]),(is,idx));
542: return(0);
543: }
545: /*@C
546: ISBlockRestoreIndices - Restores the indices associated with each block.
548: Not Collective
550: Input Parameter:
551: . is - the index set
553: Output Parameter:
554: . idx - the integer indices
556: Level: intermediate
558: Concepts: IS^block
559: Concepts: index sets^getting indices
560: Concepts: index sets^block
562: .seealso: ISRestoreIndices(), ISBlockGetIndices()
563: @*/
564: PetscErrorCode ISBlockRestoreIndices(IS is,const PetscInt *idx[])
565: {
569: PetscUseMethod(is,"ISBlockRestoreIndices_C",(IS,const PetscInt*[]),(is,idx));
570: return(0);
571: }
573: /*@
574: ISBlockGetLocalSize - Returns the local number of blocks in the index set.
576: Not Collective
578: Input Parameter:
579: . is - the index set
581: Output Parameter:
582: . size - the local number of blocks
584: Level: intermediate
586: Concepts: IS^block sizes
587: Concepts: index sets^block sizes
589: .seealso: ISGetBlockSize(), ISBlockGetSize(), ISGetSize(), ISCreateBlock()
590: @*/
591: PetscErrorCode ISBlockGetLocalSize(IS is,PetscInt *size)
592: {
596: PetscUseMethod(is,"ISBlockGetLocalSize_C",(IS,PetscInt*),(is,size));
597: return(0);
598: }
600: static PetscErrorCode ISBlockGetLocalSize_Block(IS is,PetscInt *size)
601: {
602: PetscInt bs, n;
606: PetscLayoutGetBlockSize(is->map, &bs);
607: PetscLayoutGetLocalSize(is->map, &n);
608: *size = n/bs;
609: return(0);
610: }
612: /*@
613: ISBlockGetSize - Returns the global number of blocks in the index set.
615: Not Collective
617: Input Parameter:
618: . is - the index set
620: Output Parameter:
621: . size - the global number of blocks
623: Level: intermediate
625: Concepts: IS^block sizes
626: Concepts: index sets^block sizes
628: .seealso: ISGetBlockSize(), ISBlockGetLocalSize(), ISGetSize(), ISCreateBlock()
629: @*/
630: PetscErrorCode ISBlockGetSize(IS is,PetscInt *size)
631: {
635: PetscUseMethod(is,"ISBlockGetSize_C",(IS,PetscInt*),(is,size));
636: return(0);
637: }
639: static PetscErrorCode ISBlockGetSize_Block(IS is,PetscInt *size)
640: {
641: PetscInt bs, N;
645: PetscLayoutGetBlockSize(is->map, &bs);
646: PetscLayoutGetSize(is->map, &N);
647: *size = N/bs;
648: return(0);
649: }
651: PETSC_EXTERN PetscErrorCode ISCreate_Block(IS is)
652: {
654: IS_Block *sub;
657: PetscNewLog(is,&sub);
658: is->data = (void *) sub;
659: PetscMemcpy(is->ops,&myops,sizeof(myops));
660: PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",ISBlockSetIndices_Block);
661: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",ISBlockGetIndices_Block);
662: PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",ISBlockRestoreIndices_Block);
663: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",ISBlockGetSize_Block);
664: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",ISBlockGetLocalSize_Block);
665: return(0);
666: }