Actual source code: iscoloring.c
petsc-3.10.5 2019-03-28
2: #include <petsc/private/isimpl.h>
3: #include <petscviewer.h>
4: #include <petscsf.h>
6: const char *const ISColoringTypes[] = {"global","ghosted","ISColoringType","IS_COLORING_",0};
8: PetscErrorCode ISColoringReference(ISColoring coloring)
9: {
11: (coloring)->refct++;
12: return(0);
13: }
15: PetscErrorCode ISColoringSetType(ISColoring coloring,ISColoringType type)
16: {
18: (coloring)->ctype = type;
19: return(0);
20: }
22: /*@
23: ISColoringDestroy - Destroys a coloring context.
25: Collective on ISColoring
27: Input Parameter:
28: . iscoloring - the coloring context
30: Level: advanced
32: .seealso: ISColoringView(), MatColoring
33: @*/
34: PetscErrorCode ISColoringDestroy(ISColoring *iscoloring)
35: {
36: PetscInt i;
40: if (!*iscoloring) return(0);
42: if (--(*iscoloring)->refct > 0) {*iscoloring = 0; return(0);}
44: if ((*iscoloring)->is) {
45: for (i=0; i<(*iscoloring)->n; i++) {
46: ISDestroy(&(*iscoloring)->is[i]);
47: }
48: PetscFree((*iscoloring)->is);
49: }
50: if ((*iscoloring)->allocated) {PetscFree((*iscoloring)->colors);}
51: PetscCommDestroy(&(*iscoloring)->comm);
52: PetscFree((*iscoloring));
53: return(0);
54: }
56: /*
57: ISColoringViewFromOptions - Processes command line options to determine if/how an ISColoring object is to be viewed.
59: Collective on ISColoring
61: Input Parameters:
62: + obj - the ISColoring object
63: . prefix - prefix to use for viewing, or NULL to use prefix of 'mat'
64: - optionname - option to activate viewing
66: Level: intermediate
68: Developer Note: This cannot use PetscObjectViewFromOptions() because ISColoring is not a PetscObject
70: */
71: PetscErrorCode ISColoringViewFromOptions(ISColoring obj,PetscObject bobj,const char optionname[])
72: {
73: PetscErrorCode ierr;
74: PetscViewer viewer;
75: PetscBool flg;
76: PetscViewerFormat format;
77: char *prefix;
80: prefix = bobj ? bobj->prefix : NULL;
81: PetscOptionsGetViewer(obj->comm,prefix,optionname,&viewer,&format,&flg);
82: if (flg) {
83: PetscViewerPushFormat(viewer,format);
84: ISColoringView(obj,viewer);
85: PetscViewerPopFormat(viewer);
86: PetscViewerDestroy(&viewer);
87: }
88: return(0);
89: }
91: /*@C
92: ISColoringView - Views a coloring context.
94: Collective on ISColoring
96: Input Parameters:
97: + iscoloring - the coloring context
98: - viewer - the viewer
100: Level: advanced
102: .seealso: ISColoringDestroy(), ISColoringGetIS(), MatColoring
103: @*/
104: PetscErrorCode ISColoringView(ISColoring iscoloring,PetscViewer viewer)
105: {
106: PetscInt i;
108: PetscBool iascii;
109: IS *is;
113: if (!viewer) {
114: PetscViewerASCIIGetStdout(iscoloring->comm,&viewer);
115: }
118: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
119: if (iascii) {
120: MPI_Comm comm;
121: PetscMPIInt size,rank;
123: PetscObjectGetComm((PetscObject)viewer,&comm);
124: MPI_Comm_size(comm,&size);
125: MPI_Comm_rank(comm,&rank);
126: PetscViewerASCIIPrintf(viewer,"ISColoring Object: %d MPI processes\n",size);
127: PetscViewerASCIIPushSynchronized(viewer);
128: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of colors %d\n",rank,iscoloring->n);
129: PetscViewerFlush(viewer);
130: PetscViewerASCIIPopSynchronized(viewer);
131: }
133: ISColoringGetIS(iscoloring,PETSC_IGNORE,&is);
134: for (i=0; i<iscoloring->n; i++) {
135: ISView(iscoloring->is[i],viewer);
136: }
137: ISColoringRestoreIS(iscoloring,&is);
138: return(0);
139: }
141: /*@C
142: ISColoringGetIS - Extracts index sets from the coloring context
144: Collective on ISColoring
146: Input Parameter:
147: . iscoloring - the coloring context
149: Output Parameters:
150: + nn - number of index sets in the coloring context
151: - is - array of index sets
153: Level: advanced
155: .seealso: ISColoringRestoreIS(), ISColoringView()
156: @*/
157: PetscErrorCode ISColoringGetIS(ISColoring iscoloring,PetscInt *nn,IS *isis[])
158: {
164: if (nn) *nn = iscoloring->n;
165: if (isis) {
166: if (!iscoloring->is) {
167: PetscInt *mcolors,**ii,nc = iscoloring->n,i,base, n = iscoloring->N;
168: ISColoringValue *colors = iscoloring->colors;
169: IS *is;
171: #if defined(PETSC_USE_DEBUG)
172: for (i=0; i<n; i++) {
173: if (((PetscInt)colors[i]) >= nc) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coloring is our of range index %d value %d number colors %d",(int)i,(int)colors[i],(int)nc);
174: }
175: #endif
177: /* generate the lists of nodes for each color */
178: PetscCalloc1(nc,&mcolors);
179: for (i=0; i<n; i++) mcolors[colors[i]]++;
181: PetscMalloc1(nc,&ii);
182: PetscMalloc1(n,&ii[0]);
183: for (i=1; i<nc; i++) ii[i] = ii[i-1] + mcolors[i-1];
184: PetscMemzero(mcolors,nc*sizeof(PetscInt));
186: if (iscoloring->ctype == IS_COLORING_GLOBAL) {
187: MPI_Scan(&iscoloring->N,&base,1,MPIU_INT,MPI_SUM,iscoloring->comm);
188: base -= iscoloring->N;
189: for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
190: } else if (iscoloring->ctype == IS_COLORING_LOCAL) {
191: for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i; /* local idx */
192: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this ISColoringType type");
194: PetscMalloc1(nc,&is);
195: for (i=0; i<nc; i++) {
196: ISCreateGeneral(iscoloring->comm,mcolors[i],ii[i],PETSC_COPY_VALUES,is+i);
197: }
199: iscoloring->is = is;
200: PetscFree(ii[0]);
201: PetscFree(ii);
202: PetscFree(mcolors);
203: }
204: *isis = iscoloring->is;
205: }
206: return(0);
207: }
209: /*@C
210: ISColoringRestoreIS - Restores the index sets extracted from the coloring context
212: Collective on ISColoring
214: Input Parameter:
215: + iscoloring - the coloring context
216: - is - array of index sets
218: Level: advanced
220: .seealso: ISColoringGetIS(), ISColoringView()
221: @*/
222: PetscErrorCode ISColoringRestoreIS(ISColoring iscoloring,IS *is[])
223: {
227: /* currently nothing is done here */
228: return(0);
229: }
232: /*@
233: ISColoringCreate - Generates an ISColoring context from lists (provided
234: by each processor) of colors for each node.
236: Collective on MPI_Comm
238: Input Parameters:
239: + comm - communicator for the processors creating the coloring
240: . ncolors - max color value
241: . n - number of nodes on this processor
242: . colors - array containing the colors for this processor, color numbers begin at 0.
243: - mode - see PetscCopyMode for meaning of this flag.
245: Output Parameter:
246: . iscoloring - the resulting coloring data structure
248: Options Database Key:
249: . -is_coloring_view - Activates ISColoringView()
251: Level: advanced
253: Notes:
254: By default sets coloring type to IS_COLORING_GLOBAL
256: .seealso: MatColoringCreate(), ISColoringView(), ISColoringDestroy(), ISColoringSetType()
258: @*/
259: PetscErrorCode ISColoringCreate(MPI_Comm comm,PetscInt ncolors,PetscInt n,const ISColoringValue colors[],PetscCopyMode mode,ISColoring *iscoloring)
260: {
262: PetscMPIInt size,rank,tag;
263: PetscInt base,top,i;
264: PetscInt nc,ncwork;
265: MPI_Status status;
268: if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
269: if (ncolors > 65535) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Max color value exeeds 65535 limit. This number is unrealistic. Perhaps a bug in code?\nCurrent max: %d user rewuested: %d",IS_COLORING_MAX,ncolors);
270: else SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Max color value exeeds limit. Perhaps reconfigure PETSc with --with-is-color-value-type=short?\n Current max: %d user rewuested: %d",IS_COLORING_MAX,ncolors);
271: }
272: PetscNew(iscoloring);
273: PetscCommDuplicate(comm,&(*iscoloring)->comm,&tag);
274: comm = (*iscoloring)->comm;
276: /* compute the number of the first node on my processor */
277: MPI_Comm_size(comm,&size);
279: /* should use MPI_Scan() */
280: MPI_Comm_rank(comm,&rank);
281: if (!rank) {
282: base = 0;
283: top = n;
284: } else {
285: MPI_Recv(&base,1,MPIU_INT,rank-1,tag,comm,&status);
286: top = base+n;
287: }
288: if (rank < size-1) {
289: MPI_Send(&top,1,MPIU_INT,rank+1,tag,comm);
290: }
292: /* compute the total number of colors */
293: ncwork = 0;
294: for (i=0; i<n; i++) {
295: if (ncwork < colors[i]) ncwork = colors[i];
296: }
297: ncwork++;
298: MPIU_Allreduce(&ncwork,&nc,1,MPIU_INT,MPI_MAX,comm);
299: if (nc > ncolors) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of colors passed in %D is less then the actual number of colors in array %D",ncolors,nc);
300: (*iscoloring)->n = nc;
301: (*iscoloring)->is = 0;
302: (*iscoloring)->N = n;
303: (*iscoloring)->refct = 1;
304: (*iscoloring)->ctype = IS_COLORING_GLOBAL;
305: if (mode == PETSC_COPY_VALUES) {
306: PetscMalloc1(n,&(*iscoloring)->colors);
307: PetscLogObjectMemory((PetscObject)(*iscoloring),n*sizeof(ISColoringValue));
308: PetscMemcpy((*iscoloring)->colors,colors,n*sizeof(ISColoringValue));
309: (*iscoloring)->allocated = PETSC_TRUE;
310: } else if (mode == PETSC_OWN_POINTER) {
311: (*iscoloring)->colors = (ISColoringValue*)colors;
312: (*iscoloring)->allocated = PETSC_TRUE;
313: } else {
314: (*iscoloring)->colors = (ISColoringValue*)colors;
315: (*iscoloring)->allocated = PETSC_FALSE;
316: }
317: ISColoringViewFromOptions(*iscoloring,NULL,"-is_coloring_view");
318: PetscInfo1(0,"Number of colors %D\n",nc);
319: return(0);
320: }
322: /*@
323: ISBuildTwoSided - Takes an IS that describes where we will go. Generates an IS that contains new numbers from remote or local
324: on the IS.
326: Collective on IS
328: Input Parameters
329: . ito - an IS describes where we will go. Negative target rank will be ignored
330: . toindx - an IS describes what indices should send. NULL means sending natural numbering
332: Output Parameter:
333: . rows - contains new numbers from remote or local
335: Level: advanced
337: .seealso: MatPartitioningCreate(), ISPartitioningToNumbering(), ISPartitioningCount()
339: @*/
340: PetscErrorCode ISBuildTwoSided(IS ito,IS toindx, IS *rows)
341: {
342: const PetscInt *ito_indices,*toindx_indices;
343: PetscInt *send_indices,rstart,*recv_indices,nrecvs,nsends;
344: PetscInt *tosizes,*fromsizes,i,j,*tosizes_tmp,*tooffsets_tmp,ito_ln;
345: PetscMPIInt *toranks,*fromranks,size,target_rank,*fromperm_newtoold,nto,nfrom;
346: PetscLayout isrmap;
347: MPI_Comm comm;
348: PetscSF sf;
349: PetscSFNode *iremote;
350: PetscErrorCode ierr;
353: PetscObjectGetComm((PetscObject)ito,&comm);
354: MPI_Comm_size(comm,&size);
355: ISGetLocalSize(ito,&ito_ln);
356: /* why we do not have ISGetLayout? */
357: isrmap = ito->map;
358: PetscLayoutGetRange(isrmap,&rstart,NULL);
359: ISGetIndices(ito,&ito_indices);
360: PetscCalloc2(size,&tosizes_tmp,size+1,&tooffsets_tmp);
361: for (i=0; i<ito_ln; i++) {
362: if (ito_indices[i]<0) continue;
363: #if defined(PETSC_USE_DEBUG)
364: if (ito_indices[i]>=size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"target rank %d is larger than communicator size %d ",ito_indices[i],size);
365: #endif
366: tosizes_tmp[ito_indices[i]]++;
367: }
368: nto = 0;
369: for (i=0; i<size; i++) {
370: tooffsets_tmp[i+1] = tooffsets_tmp[i]+tosizes_tmp[i];
371: if (tosizes_tmp[i]>0) nto++;
372: }
373: PetscCalloc2(nto,&toranks,2*nto,&tosizes);
374: nto = 0;
375: for (i=0; i<size; i++) {
376: if (tosizes_tmp[i]>0) {
377: toranks[nto] = i;
378: tosizes[2*nto] = tosizes_tmp[i];/* size */
379: tosizes[2*nto+1] = tooffsets_tmp[i];/* offset */
380: nto++;
381: }
382: }
383: nsends = tooffsets_tmp[size];
384: PetscCalloc1(nsends,&send_indices);
385: if (toindx) {
386: ISGetIndices(toindx,&toindx_indices);
387: }
388: for (i=0; i<ito_ln; i++) {
389: if (ito_indices[i]<0) continue;
390: target_rank = ito_indices[i];
391: send_indices[tooffsets_tmp[target_rank]] = toindx? toindx_indices[i]:(i+rstart);
392: tooffsets_tmp[target_rank]++;
393: }
394: if (toindx) {
395: ISRestoreIndices(toindx,&toindx_indices);
396: }
397: ISRestoreIndices(ito,&ito_indices);
398: PetscFree2(tosizes_tmp,tooffsets_tmp);
399: PetscCommBuildTwoSided(comm,2,MPIU_INT,nto,toranks,tosizes,&nfrom,&fromranks,&fromsizes);
400: PetscFree2(toranks,tosizes);
401: PetscCalloc1(nfrom,&fromperm_newtoold);
402: for (i=0; i<nfrom; i++) fromperm_newtoold[i] = i;
403: PetscSortMPIIntWithArray(nfrom,fromranks,fromperm_newtoold);
404: nrecvs = 0;
405: for (i=0; i<nfrom; i++) nrecvs += fromsizes[i*2];
406: PetscCalloc1(nrecvs,&recv_indices);
407: PetscCalloc1(nrecvs,&iremote);
408: nrecvs = 0;
409: for (i=0; i<nfrom; i++) {
410: for (j=0; j<fromsizes[2*fromperm_newtoold[i]]; j++) {
411: iremote[nrecvs].rank = fromranks[i];
412: iremote[nrecvs++].index = fromsizes[2*fromperm_newtoold[i]+1]+j;
413: }
414: }
415: PetscSFCreate(comm,&sf);
416: PetscSFSetGraph(sf,nsends,nrecvs,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
417: PetscSFSetType(sf,PETSCSFBASIC);
418: /* how to put a prefix ? */
419: PetscSFSetFromOptions(sf);
420: PetscSFBcastBegin(sf,MPIU_INT,send_indices,recv_indices);
421: PetscSFBcastEnd(sf,MPIU_INT,send_indices,recv_indices);
422: PetscSFDestroy(&sf);
423: PetscFree(fromranks);
424: PetscFree(fromsizes);
425: PetscFree(fromperm_newtoold);
426: PetscFree(send_indices);
427: if (rows) {
428: PetscSortInt(nrecvs,recv_indices);
429: ISCreateGeneral(comm,nrecvs,recv_indices,PETSC_OWN_POINTER,rows);
430: } else {
431: PetscFree(recv_indices);
432: }
433: return(0);
434: }
437: /*@
438: ISPartitioningToNumbering - Takes an ISPartitioning and on each processor
439: generates an IS that contains a new global node number for each index based
440: on the partitioing.
442: Collective on IS
444: Input Parameters
445: . partitioning - a partitioning as generated by MatPartitioningApply()
446: or MatPartitioningApplyND()
448: Output Parameter:
449: . is - on each processor the index set that defines the global numbers
450: (in the new numbering) for all the nodes currently (before the partitioning)
451: on that processor
453: Level: advanced
455: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningCount()
457: @*/
458: PetscErrorCode ISPartitioningToNumbering(IS part,IS *is)
459: {
460: MPI_Comm comm;
461: IS ndorder;
462: PetscInt i,np,npt,n,*starts = NULL,*sums = NULL,*lsizes = NULL,*newi = NULL;
463: const PetscInt *indices = NULL;
469: /* see if the partitioning comes from nested dissection */
470: PetscObjectQuery((PetscObject)part,"_petsc_matpartitioning_ndorder",(PetscObject*)&ndorder);
471: if (ndorder) {
472: PetscObjectReference((PetscObject)ndorder);
473: *is = ndorder;
474: return(0);
475: }
477: PetscObjectGetComm((PetscObject)part,&comm);
478: /* count the number of partitions, i.e., virtual processors */
479: ISGetLocalSize(part,&n);
480: ISGetIndices(part,&indices);
481: np = 0;
482: for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
483: MPIU_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
484: np = npt+1; /* so that it looks like a MPI_Comm_size output */
486: /*
487: lsizes - number of elements of each partition on this particular processor
488: sums - total number of "previous" nodes for any particular partition
489: starts - global number of first element in each partition on this processor
490: */
491: PetscMalloc3(np,&lsizes,np,&starts,np,&sums);
492: PetscMemzero(lsizes,np*sizeof(PetscInt));
493: for (i=0; i<n; i++) lsizes[indices[i]]++;
494: MPIU_Allreduce(lsizes,sums,np,MPIU_INT,MPI_SUM,comm);
495: MPI_Scan(lsizes,starts,np,MPIU_INT,MPI_SUM,comm);
496: for (i=0; i<np; i++) starts[i] -= lsizes[i];
497: for (i=1; i<np; i++) {
498: sums[i] += sums[i-1];
499: starts[i] += sums[i-1];
500: }
502: /*
503: For each local index give it the new global number
504: */
505: PetscMalloc1(n,&newi);
506: for (i=0; i<n; i++) newi[i] = starts[indices[i]]++;
507: PetscFree3(lsizes,starts,sums);
509: ISRestoreIndices(part,&indices);
510: ISCreateGeneral(comm,n,newi,PETSC_OWN_POINTER,is);
511: ISSetPermutation(*is);
512: return(0);
513: }
515: /*@
516: ISPartitioningCount - Takes a ISPartitioning and determines the number of
517: resulting elements on each (partition) process
519: Collective on IS
521: Input Parameters:
522: + partitioning - a partitioning as generated by MatPartitioningApply() or
523: MatPartitioningApplyND()
524: - len - length of the array count, this is the total number of partitions
526: Output Parameter:
527: . count - array of length size, to contain the number of elements assigned
528: to each partition, where size is the number of partitions generated
529: (see notes below).
531: Level: advanced
533: Notes:
534: By default the number of partitions generated (and thus the length
535: of count) is the size of the communicator associated with IS,
536: but it can be set by MatPartitioningSetNParts. The resulting array
537: of lengths can for instance serve as input of PCBJacobiSetTotalBlocks.
538: If the partitioning has been obtained by MatPartitioningApplyND(),
539: the returned count does not include the separators.
541: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningToNumbering(),
542: MatPartitioningSetNParts(), MatPartitioningApply(), MatPartitioningApplyND()
544: @*/
545: PetscErrorCode ISPartitioningCount(IS part,PetscInt len,PetscInt count[])
546: {
547: MPI_Comm comm;
548: PetscInt i,n,*lsizes;
549: const PetscInt *indices;
551: PetscMPIInt npp;
554: PetscObjectGetComm((PetscObject)part,&comm);
555: if (len == PETSC_DEFAULT) {
556: PetscMPIInt size;
557: MPI_Comm_size(comm,&size);
558: len = (PetscInt) size;
559: }
561: /* count the number of partitions */
562: ISGetLocalSize(part,&n);
563: ISGetIndices(part,&indices);
564: #if defined(PETSC_USE_DEBUG)
565: {
566: PetscInt np = 0,npt;
567: for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
568: MPIU_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
569: np = npt+1; /* so that it looks like a MPI_Comm_size output */
570: if (np > len) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Length of count array %D is less than number of partitions %D",len,np);
571: }
572: #endif
574: /*
575: lsizes - number of elements of each partition on this particular processor
576: sums - total number of "previous" nodes for any particular partition
577: starts - global number of first element in each partition on this processor
578: */
579: PetscCalloc1(len,&lsizes);
580: for (i=0; i<n; i++) {
581: if (indices[i] > -1) lsizes[indices[i]]++;
582: }
583: ISRestoreIndices(part,&indices);
584: PetscMPIIntCast(len,&npp);
585: MPIU_Allreduce(lsizes,count,npp,MPIU_INT,MPI_SUM,comm);
586: PetscFree(lsizes);
587: return(0);
588: }
590: /*@
591: ISAllGather - Given an index set (IS) on each processor, generates a large
592: index set (same on each processor) by concatenating together each
593: processors index set.
595: Collective on IS
597: Input Parameter:
598: . is - the distributed index set
600: Output Parameter:
601: . isout - the concatenated index set (same on all processors)
603: Notes:
604: ISAllGather() is clearly not scalable for large index sets.
606: The IS created on each processor must be created with a common
607: communicator (e.g., PETSC_COMM_WORLD). If the index sets were created
608: with PETSC_COMM_SELF, this routine will not work as expected, since
609: each process will generate its own new IS that consists only of
610: itself.
612: The communicator for this new IS is PETSC_COMM_SELF
614: Level: intermediate
616: Concepts: gather^index sets
617: Concepts: index sets^gathering to all processors
618: Concepts: IS^gathering to all processors
620: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock()
621: @*/
622: PetscErrorCode ISAllGather(IS is,IS *isout)
623: {
625: PetscInt *indices,n,i,N,step,first;
626: const PetscInt *lindices;
627: MPI_Comm comm;
628: PetscMPIInt size,*sizes = NULL,*offsets = NULL,nn;
629: PetscBool stride;
635: PetscObjectGetComm((PetscObject)is,&comm);
636: MPI_Comm_size(comm,&size);
637: ISGetLocalSize(is,&n);
638: PetscObjectTypeCompare((PetscObject)is,ISSTRIDE,&stride);
639: if (size == 1 && stride) { /* should handle parallel ISStride also */
640: ISStrideGetInfo(is,&first,&step);
641: ISCreateStride(PETSC_COMM_SELF,n,first,step,isout);
642: } else {
643: PetscMalloc2(size,&sizes,size,&offsets);
645: PetscMPIIntCast(n,&nn);
646: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
647: offsets[0] = 0;
648: for (i=1; i<size; i++) {
649: PetscInt s = offsets[i-1] + sizes[i-1];
650: PetscMPIIntCast(s,&offsets[i]);
651: }
652: N = offsets[size-1] + sizes[size-1];
654: PetscMalloc1(N,&indices);
655: ISGetIndices(is,&lindices);
656: MPI_Allgatherv((void*)lindices,nn,MPIU_INT,indices,sizes,offsets,MPIU_INT,comm);
657: ISRestoreIndices(is,&lindices);
658: PetscFree2(sizes,offsets);
660: ISCreateGeneral(PETSC_COMM_SELF,N,indices,PETSC_OWN_POINTER,isout);
661: }
662: return(0);
663: }
665: /*@C
666: ISAllGatherColors - Given a a set of colors on each processor, generates a large
667: set (same on each processor) by concatenating together each processors colors
669: Collective on MPI_Comm
671: Input Parameter:
672: + comm - communicator to share the indices
673: . n - local size of set
674: - lindices - local colors
676: Output Parameter:
677: + outN - total number of indices
678: - outindices - all of the colors
680: Notes:
681: ISAllGatherColors() is clearly not scalable for large index sets.
684: Level: intermediate
686: Concepts: gather^index sets
687: Concepts: index sets^gathering to all processors
688: Concepts: IS^gathering to all processors
690: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
691: @*/
692: PetscErrorCode ISAllGatherColors(MPI_Comm comm,PetscInt n,ISColoringValue *lindices,PetscInt *outN,ISColoringValue *outindices[])
693: {
694: ISColoringValue *indices;
695: PetscErrorCode ierr;
696: PetscInt i,N;
697: PetscMPIInt size,*offsets = NULL,*sizes = NULL, nn = n;
700: MPI_Comm_size(comm,&size);
701: PetscMalloc2(size,&sizes,size,&offsets);
703: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
704: offsets[0] = 0;
705: for (i=1; i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
706: N = offsets[size-1] + sizes[size-1];
707: PetscFree2(sizes,offsets);
709: PetscMalloc1(N+1,&indices);
710: MPI_Allgatherv(lindices,(PetscMPIInt)n,MPIU_COLORING_VALUE,indices,sizes,offsets,MPIU_COLORING_VALUE,comm);
712: *outindices = indices;
713: if (outN) *outN = N;
714: return(0);
715: }
717: /*@
718: ISComplement - Given an index set (IS) generates the complement index set. That is all
719: all indices that are NOT in the given set.
721: Collective on IS
723: Input Parameter:
724: + is - the index set
725: . nmin - the first index desired in the local part of the complement
726: - nmax - the largest index desired in the local part of the complement (note that all indices in is must be greater or equal to nmin and less than nmax)
728: Output Parameter:
729: . isout - the complement
731: Notes:
732: The communicator for this new IS is the same as for the input IS
734: For a parallel IS, this will generate the local part of the complement on each process
736: To generate the entire complement (on each process) of a parallel IS, first call ISAllGather() and then
737: call this routine.
739: Level: intermediate
741: Concepts: gather^index sets
742: Concepts: index sets^gathering to all processors
743: Concepts: IS^gathering to all processors
745: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
746: @*/
747: PetscErrorCode ISComplement(IS is,PetscInt nmin,PetscInt nmax,IS *isout)
748: {
750: const PetscInt *indices;
751: PetscInt n,i,j,unique,cnt,*nindices;
752: PetscBool sorted;
757: if (nmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nmin %D cannot be negative",nmin);
758: if (nmin > nmax) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nmin %D cannot be greater than nmax %D",nmin,nmax);
759: ISSorted(is,&sorted);
760: if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Index set must be sorted");
762: ISGetLocalSize(is,&n);
763: ISGetIndices(is,&indices);
764: #if defined(PETSC_USE_DEBUG)
765: for (i=0; i<n; i++) {
766: if (indices[i] < nmin) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D's value %D is smaller than minimum given %D",i,indices[i],nmin);
767: if (indices[i] >= nmax) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D's value %D is larger than maximum given %D",i,indices[i],nmax);
768: }
769: #endif
770: /* Count number of unique entries */
771: unique = (n>0);
772: for (i=0; i<n-1; i++) {
773: if (indices[i+1] != indices[i]) unique++;
774: }
775: PetscMalloc1(nmax-nmin-unique,&nindices);
776: cnt = 0;
777: for (i=nmin,j=0; i<nmax; i++) {
778: if (j<n && i==indices[j]) do { j++; } while (j<n && i==indices[j]);
779: else nindices[cnt++] = i;
780: }
781: if (cnt != nmax-nmin-unique) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number of entries found in complement %D does not match expected %D",cnt,nmax-nmin-unique);
782: ISCreateGeneral(PetscObjectComm((PetscObject)is),cnt,nindices,PETSC_OWN_POINTER,isout);
783: ISRestoreIndices(is,&indices);
784: return(0);
785: }