Actual source code: iscoloring.c
petsc-3.4.5 2014-06-29
2: #include <petsc-private/isimpl.h> /*I "petscis.h" I*/
3: #include <petscviewer.h>
5: const char *const ISColoringTypes[] = {"global","ghosted","ISColoringType","IS_COLORING_",0};
9: PetscErrorCode ISColoringReference(ISColoring coloring)
10: {
12: (coloring)->refct++;
13: return(0);
14: }
18: PetscErrorCode ISColoringSetType(ISColoring coloring,ISColoringType type)
19: {
21: (coloring)->ctype = type;
22: return(0);
23: }
27: /*@
28: ISColoringDestroy - Destroys a coloring context.
30: Collective on ISColoring
32: Input Parameter:
33: . iscoloring - the coloring context
35: Level: advanced
37: .seealso: ISColoringView(), MatGetColoring()
38: @*/
39: PetscErrorCode ISColoringDestroy(ISColoring *iscoloring)
40: {
41: PetscInt i;
45: if (!*iscoloring) return(0);
47: if (--(*iscoloring)->refct > 0) {*iscoloring = 0; return(0);}
49: if ((*iscoloring)->is) {
50: for (i=0; i<(*iscoloring)->n; i++) {
51: ISDestroy(&(*iscoloring)->is[i]);
52: }
53: PetscFree((*iscoloring)->is);
54: }
55: PetscFree((*iscoloring)->colors);
56: PetscCommDestroy(&(*iscoloring)->comm);
57: PetscFree((*iscoloring));
58: return(0);
59: }
63: /*@C
64: ISColoringView - Views a coloring context.
66: Collective on ISColoring
68: Input Parameters:
69: + iscoloring - the coloring context
70: - viewer - the viewer
72: Level: advanced
74: .seealso: ISColoringDestroy(), ISColoringGetIS(), MatGetColoring()
75: @*/
76: PetscErrorCode ISColoringView(ISColoring iscoloring,PetscViewer viewer)
77: {
78: PetscInt i;
80: PetscBool iascii;
81: IS *is;
85: if (!viewer) {
86: PetscViewerASCIIGetStdout(iscoloring->comm,&viewer);
87: }
90: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
91: if (iascii) {
92: MPI_Comm comm;
93: PetscMPIInt rank;
94: PetscObjectGetComm((PetscObject)viewer,&comm);
95: MPI_Comm_rank(comm,&rank);
96: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
97: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of colors %d\n",rank,iscoloring->n);
98: PetscViewerFlush(viewer);
99: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
100: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for ISColoring",((PetscObject)viewer)->type_name);
102: ISColoringGetIS(iscoloring,PETSC_IGNORE,&is);
103: for (i=0; i<iscoloring->n; i++) {
104: ISView(iscoloring->is[i],viewer);
105: }
106: ISColoringRestoreIS(iscoloring,&is);
107: return(0);
108: }
112: /*@C
113: ISColoringGetIS - Extracts index sets from the coloring context
115: Collective on ISColoring
117: Input Parameter:
118: . iscoloring - the coloring context
120: Output Parameters:
121: + nn - number of index sets in the coloring context
122: - is - array of index sets
124: Level: advanced
126: .seealso: ISColoringRestoreIS(), ISColoringView()
127: @*/
128: PetscErrorCode ISColoringGetIS(ISColoring iscoloring,PetscInt *nn,IS *isis[])
129: {
135: if (nn) *nn = iscoloring->n;
136: if (isis) {
137: if (!iscoloring->is) {
138: PetscInt *mcolors,**ii,nc = iscoloring->n,i,base, n = iscoloring->N;
139: ISColoringValue *colors = iscoloring->colors;
140: IS *is;
142: #if defined(PETSC_USE_DEBUG)
143: for (i=0; i<n; i++) {
144: 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);
145: }
146: #endif
148: /* generate the lists of nodes for each color */
149: PetscMalloc(nc*sizeof(PetscInt),&mcolors);
150: PetscMemzero(mcolors,nc*sizeof(PetscInt));
151: for (i=0; i<n; i++) mcolors[colors[i]]++;
153: PetscMalloc(nc*sizeof(PetscInt*),&ii);
154: PetscMalloc(n*sizeof(PetscInt),&ii[0]);
155: for (i=1; i<nc; i++) ii[i] = ii[i-1] + mcolors[i-1];
156: PetscMemzero(mcolors,nc*sizeof(PetscInt));
158: if (iscoloring->ctype == IS_COLORING_GLOBAL) {
159: MPI_Scan(&iscoloring->N,&base,1,MPIU_INT,MPI_SUM,iscoloring->comm);
160: base -= iscoloring->N;
161: for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
162: } else if (iscoloring->ctype == IS_COLORING_GHOSTED) {
163: for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i; /* local idx */
164: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this ISColoringType type");
166: PetscMalloc(nc*sizeof(IS),&is);
167: for (i=0; i<nc; i++) {
168: ISCreateGeneral(iscoloring->comm,mcolors[i],ii[i],PETSC_COPY_VALUES,is+i);
169: }
171: iscoloring->is = is;
172: PetscFree(ii[0]);
173: PetscFree(ii);
174: PetscFree(mcolors);
175: }
176: *isis = iscoloring->is;
177: }
178: return(0);
179: }
183: /*@C
184: ISColoringRestoreIS - Restores the index sets extracted from the coloring context
186: Collective on ISColoring
188: Input Parameter:
189: + iscoloring - the coloring context
190: - is - array of index sets
192: Level: advanced
194: .seealso: ISColoringGetIS(), ISColoringView()
195: @*/
196: PetscErrorCode ISColoringRestoreIS(ISColoring iscoloring,IS *is[])
197: {
201: /* currently nothing is done here */
202: return(0);
203: }
208: /*@C
209: ISColoringCreate - Generates an ISColoring context from lists (provided
210: by each processor) of colors for each node.
212: Collective on MPI_Comm
214: Input Parameters:
215: + comm - communicator for the processors creating the coloring
216: . ncolors - max color value
217: . n - number of nodes on this processor
218: - colors - array containing the colors for this processor, color
219: numbers begin at 0. In C/C++ this array must have been obtained with PetscMalloc()
220: and should NOT be freed (The ISColoringDestroy() will free it).
222: Output Parameter:
223: . iscoloring - the resulting coloring data structure
225: Options Database Key:
226: . -is_coloring_view - Activates ISColoringView()
228: Level: advanced
230: Notes: By default sets coloring type to IS_COLORING_GLOBAL
232: .seealso: MatColoringCreate(), ISColoringView(), ISColoringDestroy(), ISColoringSetType()
234: @*/
235: PetscErrorCode ISColoringCreate(MPI_Comm comm,PetscInt ncolors,PetscInt n,const ISColoringValue colors[],ISColoring *iscoloring)
236: {
238: PetscMPIInt size,rank,tag;
239: PetscInt base,top,i;
240: PetscInt nc,ncwork;
241: PetscBool flg = PETSC_FALSE;
242: MPI_Status status;
245: if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
246: 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);
247: 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);
248: }
249: PetscNew(struct _n_ISColoring,iscoloring);
250: PetscCommDuplicate(comm,&(*iscoloring)->comm,&tag);
251: comm = (*iscoloring)->comm;
253: /* compute the number of the first node on my processor */
254: MPI_Comm_size(comm,&size);
256: /* should use MPI_Scan() */
257: MPI_Comm_rank(comm,&rank);
258: if (!rank) {
259: base = 0;
260: top = n;
261: } else {
262: MPI_Recv(&base,1,MPIU_INT,rank-1,tag,comm,&status);
263: top = base+n;
264: }
265: if (rank < size-1) {
266: MPI_Send(&top,1,MPIU_INT,rank+1,tag,comm);
267: }
269: /* compute the total number of colors */
270: ncwork = 0;
271: for (i=0; i<n; i++) {
272: if (ncwork < colors[i]) ncwork = colors[i];
273: }
274: ncwork++;
275: MPI_Allreduce(&ncwork,&nc,1,MPIU_INT,MPI_MAX,comm);
276: 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);
277: (*iscoloring)->n = nc;
278: (*iscoloring)->is = 0;
279: (*iscoloring)->colors = (ISColoringValue*)colors;
280: (*iscoloring)->N = n;
281: (*iscoloring)->refct = 1;
282: (*iscoloring)->ctype = IS_COLORING_GLOBAL;
284: PetscOptionsGetBool(NULL,"-is_coloring_view",&flg,NULL);
285: if (flg) {
286: PetscViewer viewer;
287: PetscViewerASCIIGetStdout((*iscoloring)->comm,&viewer);
288: ISColoringView(*iscoloring,viewer);
289: }
290: PetscInfo1(0,"Number of colors %D\n",nc);
291: return(0);
292: }
296: /*@
297: ISPartitioningToNumbering - Takes an ISPartitioning and on each processor
298: generates an IS that contains a new global node number for each index based
299: on the partitioing.
301: Collective on IS
303: Input Parameters
304: . partitioning - a partitioning as generated by MatPartitioningApply()
306: Output Parameter:
307: . is - on each processor the index set that defines the global numbers
308: (in the new numbering) for all the nodes currently (before the partitioning)
309: on that processor
311: Level: advanced
313: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningCount()
315: @*/
316: PetscErrorCode ISPartitioningToNumbering(IS part,IS *is)
317: {
318: MPI_Comm comm;
319: PetscInt i,np,npt,n,*starts = NULL,*sums = NULL,*lsizes = NULL,*newi = NULL;
320: const PetscInt *indices = NULL;
324: PetscObjectGetComm((PetscObject)part,&comm);
326: /* count the number of partitions, i.e., virtual processors */
327: ISGetLocalSize(part,&n);
328: ISGetIndices(part,&indices);
329: np = 0;
330: for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
331: MPI_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
332: np = npt+1; /* so that it looks like a MPI_Comm_size output */
334: /*
335: lsizes - number of elements of each partition on this particular processor
336: sums - total number of "previous" nodes for any particular partition
337: starts - global number of first element in each partition on this processor
338: */
339: PetscMalloc3(np,PetscInt,&lsizes,np,PetscInt,&starts,np,PetscInt,&sums);
340: PetscMemzero(lsizes,np*sizeof(PetscInt));
341: for (i=0; i<n; i++) lsizes[indices[i]]++;
342: MPI_Allreduce(lsizes,sums,np,MPIU_INT,MPI_SUM,comm);
343: MPI_Scan(lsizes,starts,np,MPIU_INT,MPI_SUM,comm);
344: for (i=0; i<np; i++) starts[i] -= lsizes[i];
345: for (i=1; i<np; i++) {
346: sums[i] += sums[i-1];
347: starts[i] += sums[i-1];
348: }
350: /*
351: For each local index give it the new global number
352: */
353: PetscMalloc(n*sizeof(PetscInt),&newi);
354: for (i=0; i<n; i++) newi[i] = starts[indices[i]]++;
355: PetscFree3(lsizes,starts,sums);
357: ISRestoreIndices(part,&indices);
358: ISCreateGeneral(comm,n,newi,PETSC_OWN_POINTER,is);
359: ISSetPermutation(*is);
360: return(0);
361: }
365: /*@
366: ISPartitioningCount - Takes a ISPartitioning and determines the number of
367: resulting elements on each (partition) process
369: Collective on IS
371: Input Parameters:
372: + partitioning - a partitioning as generated by MatPartitioningApply()
373: - len - length of the array count, this is the total number of partitions
375: Output Parameter:
376: . count - array of length size, to contain the number of elements assigned
377: to each partition, where size is the number of partitions generated
378: (see notes below).
380: Level: advanced
382: Notes:
383: By default the number of partitions generated (and thus the length
384: of count) is the size of the communicator associated with IS,
385: but it can be set by MatPartitioningSetNParts. The resulting array
386: of lengths can for instance serve as input of PCBJacobiSetTotalBlocks.
389: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningToNumbering(),
390: MatPartitioningSetNParts(), MatPartitioningApply()
392: @*/
393: PetscErrorCode ISPartitioningCount(IS part,PetscInt len,PetscInt count[])
394: {
395: MPI_Comm comm;
396: PetscInt i,n,*lsizes;
397: const PetscInt *indices;
399: PetscMPIInt npp;
402: PetscObjectGetComm((PetscObject)part,&comm);
403: if (len == PETSC_DEFAULT) {
404: PetscMPIInt size;
405: MPI_Comm_size(comm,&size);
406: len = (PetscInt) size;
407: }
409: /* count the number of partitions */
410: ISGetLocalSize(part,&n);
411: ISGetIndices(part,&indices);
412: #if defined(PETSC_USE_DEBUG)
413: {
414: PetscInt np = 0,npt;
415: for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
416: MPI_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
417: np = npt+1; /* so that it looks like a MPI_Comm_size output */
418: 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);
419: }
420: #endif
422: /*
423: lsizes - number of elements of each partition on this particular processor
424: sums - total number of "previous" nodes for any particular partition
425: starts - global number of first element in each partition on this processor
426: */
427: PetscMalloc(len*sizeof(PetscInt),&lsizes);
428: PetscMemzero(lsizes,len*sizeof(PetscInt));
429: for (i=0; i<n; i++) lsizes[indices[i]]++;
430: ISRestoreIndices(part,&indices);
431: PetscMPIIntCast(len,&npp);
432: MPI_Allreduce(lsizes,count,npp,MPIU_INT,MPI_SUM,comm);
433: PetscFree(lsizes);
434: return(0);
435: }
439: /*@
440: ISAllGather - Given an index set (IS) on each processor, generates a large
441: index set (same on each processor) by concatenating together each
442: processors index set.
444: Collective on IS
446: Input Parameter:
447: . is - the distributed index set
449: Output Parameter:
450: . isout - the concatenated index set (same on all processors)
452: Notes:
453: ISAllGather() is clearly not scalable for large index sets.
455: The IS created on each processor must be created with a common
456: communicator (e.g., PETSC_COMM_WORLD). If the index sets were created
457: with PETSC_COMM_SELF, this routine will not work as expected, since
458: each process will generate its own new IS that consists only of
459: itself.
461: The communicator for this new IS is PETSC_COMM_SELF
463: Level: intermediate
465: Concepts: gather^index sets
466: Concepts: index sets^gathering to all processors
467: Concepts: IS^gathering to all processors
469: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock()
470: @*/
471: PetscErrorCode ISAllGather(IS is,IS *isout)
472: {
474: PetscInt *indices,n,i,N,step,first;
475: const PetscInt *lindices;
476: MPI_Comm comm;
477: PetscMPIInt size,*sizes = NULL,*offsets = NULL,nn;
478: PetscBool stride;
484: PetscObjectGetComm((PetscObject)is,&comm);
485: MPI_Comm_size(comm,&size);
486: ISGetLocalSize(is,&n);
487: PetscObjectTypeCompare((PetscObject)is,ISSTRIDE,&stride);
488: if (size == 1 && stride) { /* should handle parallel ISStride also */
489: ISStrideGetInfo(is,&first,&step);
490: ISCreateStride(PETSC_COMM_SELF,n,first,step,isout);
491: } else {
492: PetscMalloc2(size,PetscMPIInt,&sizes,size,PetscMPIInt,&offsets);
494: PetscMPIIntCast(n,&nn);
495: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
496: offsets[0] = 0;
497: for (i=1; i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
498: N = offsets[size-1] + sizes[size-1];
500: PetscMalloc(N*sizeof(PetscInt),&indices);
501: ISGetIndices(is,&lindices);
502: MPI_Allgatherv((void*)lindices,nn,MPIU_INT,indices,sizes,offsets,MPIU_INT,comm);
503: ISRestoreIndices(is,&lindices);
504: PetscFree2(sizes,offsets);
506: ISCreateGeneral(PETSC_COMM_SELF,N,indices,PETSC_OWN_POINTER,isout);
507: }
508: return(0);
509: }
513: /*@C
514: ISAllGatherColors - Given a a set of colors on each processor, generates a large
515: set (same on each processor) by concatenating together each processors colors
517: Collective on MPI_Comm
519: Input Parameter:
520: + comm - communicator to share the indices
521: . n - local size of set
522: - lindices - local colors
524: Output Parameter:
525: + outN - total number of indices
526: - outindices - all of the colors
528: Notes:
529: ISAllGatherColors() is clearly not scalable for large index sets.
532: Level: intermediate
534: Concepts: gather^index sets
535: Concepts: index sets^gathering to all processors
536: Concepts: IS^gathering to all processors
538: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
539: @*/
540: PetscErrorCode ISAllGatherColors(MPI_Comm comm,PetscInt n,ISColoringValue *lindices,PetscInt *outN,ISColoringValue *outindices[])
541: {
542: ISColoringValue *indices;
543: PetscErrorCode ierr;
544: PetscInt i,N;
545: PetscMPIInt size,*offsets = NULL,*sizes = NULL, nn = n;
548: MPI_Comm_size(comm,&size);
549: PetscMalloc2(size,PetscMPIInt,&sizes,size,PetscMPIInt,&offsets);
551: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
552: offsets[0] = 0;
553: for (i=1; i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
554: N = offsets[size-1] + sizes[size-1];
555: PetscFree2(sizes,offsets);
557: PetscMalloc((N+1)*sizeof(ISColoringValue),&indices);
558: MPI_Allgatherv(lindices,(PetscMPIInt)n,MPIU_COLORING_VALUE,indices,sizes,offsets,MPIU_COLORING_VALUE,comm);
560: *outindices = indices;
561: if (outN) *outN = N;
562: return(0);
563: }
567: /*@
568: ISComplement - Given an index set (IS) generates the complement index set. That is all
569: all indices that are NOT in the given set.
571: Collective on IS
573: Input Parameter:
574: + is - the index set
575: . nmin - the first index desired in the local part of the complement
576: - 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)
578: Output Parameter:
579: . isout - the complement
581: Notes: The communicator for this new IS is the same as for the input IS
583: For a parallel IS, this will generate the local part of the complement on each process
585: To generate the entire complement (on each process) of a parallel IS, first call ISAllGather() and then
586: call this routine.
588: Level: intermediate
590: Concepts: gather^index sets
591: Concepts: index sets^gathering to all processors
592: Concepts: IS^gathering to all processors
594: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
595: @*/
596: PetscErrorCode ISComplement(IS is,PetscInt nmin,PetscInt nmax,IS *isout)
597: {
599: const PetscInt *indices;
600: PetscInt n,i,j,unique,cnt,*nindices;
601: PetscBool sorted;
606: if (nmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nmin %D cannot be negative",nmin);
607: if (nmin > nmax) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nmin %D cannot be greater than nmax %D",nmin,nmax);
608: ISSorted(is,&sorted);
609: if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Index set must be sorted");
611: ISGetLocalSize(is,&n);
612: ISGetIndices(is,&indices);
613: #if defined(PETSC_USE_DEBUG)
614: for (i=0; i<n; i++) {
615: 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);
616: 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);
617: }
618: #endif
619: /* Count number of unique entries */
620: unique = (n>0);
621: for (i=0; i<n-1; i++) {
622: if (indices[i+1] != indices[i]) unique++;
623: }
624: PetscMalloc((nmax-nmin-unique)*sizeof(PetscInt),&nindices);
625: cnt = 0;
626: for (i=nmin,j=0; i<nmax; i++) {
627: if (j<n && i==indices[j]) do { j++; } while (j<n && i==indices[j]);
628: else nindices[cnt++] = i;
629: }
630: 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);
631: ISCreateGeneral(PetscObjectComm((PetscObject)is),cnt,nindices,PETSC_OWN_POINTER,isout);
632: ISRestoreIndices(is,&indices);
633: return(0);
634: }