Actual source code: block.c
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 <petscviewer.h>
9: typedef struct {
10: PetscBool sorted; /* are the blocks sorted? */
11: PetscBool allocated; /* did we allocate the index array ourselves? */
12: PetscInt *idx;
13: } IS_Block;
15: static PetscErrorCode ISDestroy_Block(IS is)
16: {
17: IS_Block *sub = (IS_Block*)is->data;
19: if (sub->allocated) PetscFree(sub->idx);
20: PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",NULL);
21: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",NULL);
22: PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",NULL);
23: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",NULL);
24: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",NULL);
25: PetscFree(is->data);
26: return 0;
27: }
29: static PetscErrorCode ISLocate_Block(IS is,PetscInt key,PetscInt *location)
30: {
31: IS_Block *sub = (IS_Block*)is->data;
32: PetscInt numIdx, i, bs, bkey, mkey;
33: PetscBool sorted;
35: PetscLayoutGetBlockSize(is->map,&bs);
36: PetscLayoutGetSize(is->map,&numIdx);
37: numIdx /= bs;
38: bkey = key / bs;
39: mkey = key % bs;
40: if (mkey < 0) {
41: bkey--;
42: mkey += bs;
43: }
44: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sorted);
45: if (sorted) {
46: PetscFindInt(bkey,numIdx,sub->idx,location);
47: } else {
48: const PetscInt *idx = sub->idx;
50: *location = -1;
51: for (i = 0; i < numIdx; i++) {
52: if (idx[i] == bkey) {
53: *location = i;
54: break;
55: }
56: }
57: }
58: if (*location >= 0) {
59: *location = *location * bs + mkey;
60: }
61: return 0;
62: }
64: static PetscErrorCode ISGetIndices_Block(IS in,const PetscInt *idx[])
65: {
66: IS_Block *sub = (IS_Block*)in->data;
67: PetscInt i,j,k,bs,n,*ii,*jj;
69: PetscLayoutGetBlockSize(in->map, &bs);
70: PetscLayoutGetLocalSize(in->map, &n);
71: n /= bs;
72: if (bs == 1) *idx = sub->idx;
73: else {
74: if (n) {
75: PetscMalloc1(bs*n,&jj);
76: *idx = jj;
77: k = 0;
78: ii = sub->idx;
79: for (i=0; i<n; i++)
80: for (j=0; j<bs; j++)
81: jj[k++] = bs*ii[i] + j;
82: } else {
83: /* do not malloc for zero size because F90Array1dCreate() inside ISRestoreArrayF90() does not keep array when zero length array */
84: *idx = NULL;
85: }
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;
95: PetscLayoutGetBlockSize(is->map, &bs);
96: if (bs != 1) {
97: PetscFree(*(void**)idx);
98: } else {
99: /* F90Array1dCreate() inside ISRestoreArrayF90() does not keep array when zero length array */
101: }
102: return 0;
103: }
105: static PetscErrorCode ISInvertPermutation_Block(IS is,PetscInt nlocal,IS *isout)
106: {
107: IS_Block *sub = (IS_Block*)is->data;
108: PetscInt i,*ii,bs,n,*idx = sub->idx;
109: 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: }
124: static PetscErrorCode ISView_Block(IS is, PetscViewer viewer)
125: {
126: IS_Block *sub = (IS_Block*)is->data;
127: PetscInt i,bs,n,*idx = sub->idx;
128: PetscBool iascii,ibinary;
130: PetscLayoutGetBlockSize(is->map, &bs);
131: PetscLayoutGetLocalSize(is->map, &n);
132: n /= bs;
133: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
134: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&ibinary);
135: if (iascii) {
136: PetscViewerFormat fmt;
138: PetscViewerGetFormat(viewer,&fmt);
139: if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
140: IS ist;
141: const char *name;
142: const PetscInt *idx;
143: PetscInt n;
145: PetscObjectGetName((PetscObject)is,&name);
146: ISGetLocalSize(is,&n);
147: ISGetIndices(is,&idx);
148: ISCreateGeneral(PetscObjectComm((PetscObject)is),n,idx,PETSC_USE_POINTER,&ist);
149: PetscObjectSetName((PetscObject)ist,name);
150: ISView(ist,viewer);
151: ISDestroy(&ist);
152: ISRestoreIndices(is,&idx);
153: } else {
154: PetscBool isperm;
156: ISGetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_FALSE,&isperm);
157: if (isperm) PetscViewerASCIIPrintf(viewer,"Block Index set is permutation\n");
158: PetscViewerASCIIPushSynchronized(viewer);
159: PetscViewerASCIISynchronizedPrintf(viewer,"Block size %" PetscInt_FMT "\n",bs);
160: PetscViewerASCIISynchronizedPrintf(viewer,"Number of block indices in set %" PetscInt_FMT "\n",n);
161: PetscViewerASCIISynchronizedPrintf(viewer,"The first indices of each block are\n");
162: for (i=0; i<n; i++) {
163: PetscViewerASCIISynchronizedPrintf(viewer,"Block %" PetscInt_FMT " Index %" PetscInt_FMT "\n",i,idx[i]);
164: }
165: PetscViewerFlush(viewer);
166: PetscViewerASCIIPopSynchronized(viewer);
167: }
168: } else if (ibinary) {
169: ISView_Binary(is,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;
179: PetscLayoutGetBlockSize(is->map, &bs);
180: PetscLayoutGetLocalSize(is->map, &n);
181: PetscIntSortSemiOrdered(n/bs,sub->idx);
182: return 0;
183: }
185: static PetscErrorCode ISSortRemoveDups_Block(IS is)
186: {
187: IS_Block *sub = (IS_Block*)is->data;
188: PetscInt bs, n, nb;
189: PetscBool sorted;
191: PetscLayoutGetBlockSize(is->map, &bs);
192: PetscLayoutGetLocalSize(is->map, &n);
193: nb = n/bs;
194: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sorted);
195: if (sorted) {
196: PetscSortedRemoveDupsInt(&nb,sub->idx);
197: } else {
198: PetscSortRemoveDupsInt(&nb,sub->idx);
199: }
200: PetscLayoutDestroy(&is->map);
201: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), nb*bs, PETSC_DECIDE, bs, &is->map);
202: return 0;
203: }
205: static PetscErrorCode ISSorted_Block(IS is,PetscBool *flg)
206: {
207: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,flg);
208: return 0;
209: }
211: static PetscErrorCode ISSortedLocal_Block(IS is,PetscBool *flg)
212: {
213: IS_Block *sub = (IS_Block*)is->data;
214: PetscInt n, bs, i, *idx;
216: PetscLayoutGetLocalSize(is->map, &n);
217: PetscLayoutGetBlockSize(is->map, &bs);
218: n /= bs;
219: idx = sub->idx;
220: for (i = 1; i < n; i++) if (idx[i] < idx[i - 1]) break;
221: if (i < n) *flg = PETSC_FALSE;
222: else *flg = PETSC_TRUE;
223: return 0;
224: }
226: static PetscErrorCode ISUniqueLocal_Block(IS is,PetscBool *flg)
227: {
228: IS_Block *sub = (IS_Block*)is->data;
229: PetscInt n, bs, i, *idx, *idxcopy = NULL;
230: PetscBool sortedLocal;
232: PetscLayoutGetLocalSize(is->map, &n);
233: PetscLayoutGetBlockSize(is->map, &bs);
234: n /= bs;
235: idx = sub->idx;
236: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sortedLocal);
237: if (!sortedLocal) {
238: PetscMalloc1(n, &idxcopy);
239: PetscArraycpy(idxcopy, idx, n);
240: PetscIntSortSemiOrdered(n, idxcopy);
241: idx = idxcopy;
242: }
243: for (i = 1; i < n; i++) if (idx[i] == idx[i - 1]) break;
244: if (i < n) *flg = PETSC_FALSE;
245: else *flg = PETSC_TRUE;
246: PetscFree(idxcopy);
247: return 0;
248: }
250: static PetscErrorCode ISPermutationLocal_Block(IS is,PetscBool *flg)
251: {
252: IS_Block *sub = (IS_Block*)is->data;
253: PetscInt n, bs, i, *idx, *idxcopy = NULL;
254: PetscBool sortedLocal;
256: PetscLayoutGetLocalSize(is->map, &n);
257: PetscLayoutGetBlockSize(is->map, &bs);
258: n /= bs;
259: idx = sub->idx;
260: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sortedLocal);
261: if (!sortedLocal) {
262: PetscMalloc1(n, &idxcopy);
263: PetscArraycpy(idxcopy, idx, n);
264: PetscIntSortSemiOrdered(n, idxcopy);
265: idx = idxcopy;
266: }
267: for (i = 0; i < n; i++) if (idx[i] != i) break;
268: if (i < n) *flg = PETSC_FALSE;
269: else *flg = PETSC_TRUE;
270: PetscFree(idxcopy);
271: return 0;
272: }
274: static PetscErrorCode ISIntervalLocal_Block(IS is,PetscBool *flg)
275: {
276: IS_Block *sub = (IS_Block*)is->data;
277: PetscInt n, bs, i, *idx;
279: PetscLayoutGetLocalSize(is->map, &n);
280: PetscLayoutGetBlockSize(is->map, &bs);
281: n /= bs;
282: idx = sub->idx;
283: for (i = 1; i < n; i++) if (idx[i] != idx[i - 1] + 1) break;
284: if (i < n) *flg = PETSC_FALSE;
285: else *flg = PETSC_TRUE;
286: return 0;
287: }
289: static PetscErrorCode ISDuplicate_Block(IS is,IS *newIS)
290: {
291: IS_Block *sub = (IS_Block*)is->data;
292: PetscInt bs, n;
294: PetscLayoutGetBlockSize(is->map, &bs);
295: PetscLayoutGetLocalSize(is->map, &n);
296: n /= bs;
297: ISCreateBlock(PetscObjectComm((PetscObject)is),bs,n,sub->idx,PETSC_COPY_VALUES,newIS);
298: return 0;
299: }
301: static PetscErrorCode ISCopy_Block(IS is,IS isy)
302: {
303: IS_Block *is_block = (IS_Block*)is->data,*isy_block = (IS_Block*)isy->data;
304: PetscInt bs, n, N, bsy, ny, Ny;
306: PetscLayoutGetBlockSize(is->map, &bs);
307: PetscLayoutGetLocalSize(is->map, &n);
308: PetscLayoutGetSize(is->map, &N);
309: PetscLayoutGetBlockSize(isy->map, &bsy);
310: PetscLayoutGetLocalSize(isy->map, &ny);
311: PetscLayoutGetSize(isy->map, &Ny);
313: PetscArraycpy(isy_block->idx,is_block->idx,(n/bs));
314: return 0;
315: }
317: static PetscErrorCode ISOnComm_Block(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
318: {
319: IS_Block *sub = (IS_Block*)is->data;
320: PetscInt bs, n;
323: PetscLayoutGetBlockSize(is->map, &bs);
324: PetscLayoutGetLocalSize(is->map, &n);
325: ISCreateBlock(comm,bs,n/bs,sub->idx,mode,newis);
326: return 0;
327: }
329: static PetscErrorCode ISSetBlockSize_Block(IS is,PetscInt bs)
330: {
332: PetscLayoutSetBlockSize(is->map, bs);
333: return 0;
334: }
336: static PetscErrorCode ISToGeneral_Block(IS inis)
337: {
338: IS_Block *sub = (IS_Block*)inis->data;
339: PetscInt bs,n;
340: const PetscInt *idx;
342: ISGetBlockSize(inis,&bs);
343: ISGetLocalSize(inis,&n);
344: ISGetIndices(inis,&idx);
345: if (bs == 1) {
346: PetscCopyMode mode = sub->allocated ? PETSC_OWN_POINTER : PETSC_USE_POINTER;
347: sub->allocated = PETSC_FALSE; /* prevent deallocation when changing the subtype*/
348: ISSetType(inis,ISGENERAL);
349: ISGeneralSetIndices(inis,n,idx,mode);
350: } else {
351: ISSetType(inis,ISGENERAL);
352: ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
353: }
354: return 0;
355: }
357: static struct _ISOps myops = { ISGetIndices_Block,
358: ISRestoreIndices_Block,
359: ISInvertPermutation_Block,
360: ISSort_Block,
361: ISSortRemoveDups_Block,
362: ISSorted_Block,
363: ISDuplicate_Block,
364: ISDestroy_Block,
365: ISView_Block,
366: ISLoad_Default,
367: ISCopy_Block,
368: ISToGeneral_Block,
369: ISOnComm_Block,
370: ISSetBlockSize_Block,
371: NULL,
372: ISLocate_Block,
373: /* we can have specialized local routines for determining properties,
374: * but unless the block size is the same on each process (which is not guaranteed at
375: * the moment), then trying to do something specialized for global properties is too
376: * complicated */
377: ISSortedLocal_Block,
378: NULL,
379: ISUniqueLocal_Block,
380: NULL,
381: ISPermutationLocal_Block,
382: NULL,
383: ISIntervalLocal_Block,
384: NULL};
386: /*@
387: ISBlockSetIndices - Set integers representing blocks of indices in an index set.
389: Collective on IS
391: Input Parameters:
392: + is - the index set
393: . bs - number of elements in each block
394: . n - the length of the index set (the number of blocks)
395: . idx - the list of integers, one for each block, the integers contain the index of the first index of each block divided by the block size
396: - mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported
398: Notes:
399: When the communicator is not MPI_COMM_SELF, the operations on the
400: index sets, IS, are NOT conceptually the same as MPI_Group operations.
401: The index sets are then distributed sets of indices and thus certain operations
402: on them are collective.
404: Example:
405: If you wish to index the values {0,1,4,5}, then use
406: a block size of 2 and idx of {0,2}.
408: Level: beginner
410: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather(), ISCreateBlock(), ISBLOCK, ISGeneralSetIndices()
411: @*/
412: PetscErrorCode ISBlockSetIndices(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
413: {
414: ISClearInfoCache(is,PETSC_FALSE);
415: PetscUseMethod(is,"ISBlockSetIndices_C",(IS,PetscInt,PetscInt,const PetscInt[],PetscCopyMode),(is,bs,n,idx,mode));
416: return 0;
417: }
419: static PetscErrorCode ISBlockSetIndices_Block(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
420: {
421: PetscInt i,min,max;
422: IS_Block *sub = (IS_Block*)is->data;
423: PetscLayout map;
429: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is),n*bs,is->map->N,bs,&map);
430: PetscLayoutDestroy(&is->map);
431: is->map = map;
433: if (sub->allocated) PetscFree(sub->idx);
434: if (mode == PETSC_COPY_VALUES) {
435: PetscMalloc1(n,&sub->idx);
436: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
437: PetscArraycpy(sub->idx,idx,n);
438: sub->allocated = PETSC_TRUE;
439: } else if (mode == PETSC_OWN_POINTER) {
440: sub->idx = (PetscInt*) idx;
441: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
442: sub->allocated = PETSC_TRUE;
443: } else if (mode == PETSC_USE_POINTER) {
444: sub->idx = (PetscInt*) idx;
445: sub->allocated = PETSC_FALSE;
446: }
448: if (n) {
449: min = max = idx[0];
450: for (i=1; i<n; i++) {
451: if (idx[i] < min) min = idx[i];
452: if (idx[i] > max) max = idx[i];
453: }
454: is->min = bs*min;
455: is->max = bs*max+bs-1;
456: } else {
457: is->min = PETSC_MAX_INT;
458: is->max = PETSC_MIN_INT;
459: }
460: return 0;
461: }
463: /*@
464: ISCreateBlock - Creates a data structure for an index set containing
465: a list of integers. Each integer represents a fixed block size set of indices.
467: Collective
469: Input Parameters:
470: + comm - the MPI communicator
471: . bs - number of elements in each block
472: . n - the length of the index set (the number of blocks)
473: . idx - the list of integers, one for each block, the integers contain the index of the first entry of each block divided by the block size
474: - mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported in this routine
476: Output Parameter:
477: . is - the new index set
479: Notes:
480: When the communicator is not MPI_COMM_SELF, the operations on the
481: index sets, IS, are NOT conceptually the same as MPI_Group operations.
482: The index sets are then distributed sets of indices and thus certain operations
483: on them are collective.
485: Example:
486: If you wish to index the values {0,1,6,7}, then use
487: a block size of 2 and idx of {0,3}.
489: Level: beginner
491: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather(), ISBlockSetIndices(), ISBLOCK, ISGENERAL
492: @*/
493: PetscErrorCode ISCreateBlock(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
494: {
500: ISCreate(comm,is);
501: ISSetType(*is,ISBLOCK);
502: ISBlockSetIndices(*is,bs,n,idx,mode);
503: return 0;
504: }
506: static PetscErrorCode ISBlockGetIndices_Block(IS is,const PetscInt *idx[])
507: {
508: IS_Block *sub = (IS_Block*)is->data;
510: *idx = sub->idx;
511: return 0;
512: }
514: static PetscErrorCode ISBlockRestoreIndices_Block(IS is,const PetscInt *idx[])
515: {
516: return 0;
517: }
519: /*@C
520: ISBlockGetIndices - Gets the indices associated with each block.
522: Not Collective
524: Input Parameter:
525: . is - the index set
527: Output Parameter:
528: . idx - the integer indices, one for each block and count of block not indices
530: Level: intermediate
532: .seealso: ISGetIndices(), ISBlockRestoreIndices(), ISBLOCK, ISBlockSetIndices(), ISCreateBlock()
533: @*/
534: PetscErrorCode ISBlockGetIndices(IS is,const PetscInt *idx[])
535: {
536: PetscUseMethod(is,"ISBlockGetIndices_C",(IS,const PetscInt*[]),(is,idx));
537: return 0;
538: }
540: /*@C
541: ISBlockRestoreIndices - Restores the indices associated with each block.
543: Not Collective
545: Input Parameter:
546: . is - the index set
548: Output Parameter:
549: . idx - the integer indices
551: Level: intermediate
553: .seealso: ISRestoreIndices(), ISBlockGetIndices()
554: @*/
555: PetscErrorCode ISBlockRestoreIndices(IS is,const PetscInt *idx[])
556: {
557: PetscUseMethod(is,"ISBlockRestoreIndices_C",(IS,const PetscInt*[]),(is,idx));
558: return 0;
559: }
561: /*@
562: ISBlockGetLocalSize - Returns the local number of blocks in the index set.
564: Not Collective
566: Input Parameter:
567: . is - the index set
569: Output Parameter:
570: . size - the local number of blocks
572: Level: intermediate
574: .seealso: ISGetBlockSize(), ISBlockGetSize(), ISGetSize(), ISCreateBlock(), ISBLOCK
575: @*/
576: PetscErrorCode ISBlockGetLocalSize(IS is,PetscInt *size)
577: {
578: PetscUseMethod(is,"ISBlockGetLocalSize_C",(IS,PetscInt*),(is,size));
579: return 0;
580: }
582: static PetscErrorCode ISBlockGetLocalSize_Block(IS is,PetscInt *size)
583: {
584: PetscInt bs, n;
586: PetscLayoutGetBlockSize(is->map, &bs);
587: PetscLayoutGetLocalSize(is->map, &n);
588: *size = n/bs;
589: return 0;
590: }
592: /*@
593: ISBlockGetSize - Returns the global number of blocks in the index set.
595: Not Collective
597: Input Parameter:
598: . is - the index set
600: Output Parameter:
601: . size - the global number of blocks
603: Level: intermediate
605: .seealso: ISGetBlockSize(), ISBlockGetLocalSize(), ISGetSize(), ISCreateBlock(), ISBLOCK
606: @*/
607: PetscErrorCode ISBlockGetSize(IS is,PetscInt *size)
608: {
609: PetscUseMethod(is,"ISBlockGetSize_C",(IS,PetscInt*),(is,size));
610: return 0;
611: }
613: static PetscErrorCode ISBlockGetSize_Block(IS is,PetscInt *size)
614: {
615: PetscInt bs, N;
617: PetscLayoutGetBlockSize(is->map, &bs);
618: PetscLayoutGetSize(is->map, &N);
619: *size = N/bs;
620: return 0;
621: }
623: PETSC_EXTERN PetscErrorCode ISCreate_Block(IS is)
624: {
625: IS_Block *sub;
627: PetscNewLog(is,&sub);
628: is->data = (void *) sub;
629: PetscMemcpy(is->ops,&myops,sizeof(myops));
630: PetscObjectComposeFunction((PetscObject)is,"ISBlockSetIndices_C",ISBlockSetIndices_Block);
631: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetIndices_C",ISBlockGetIndices_Block);
632: PetscObjectComposeFunction((PetscObject)is,"ISBlockRestoreIndices_C",ISBlockRestoreIndices_Block);
633: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetSize_C",ISBlockGetSize_Block);
634: PetscObjectComposeFunction((PetscObject)is,"ISBlockGetLocalSize_C",ISBlockGetLocalSize_Block);
635: return 0;
636: }