Actual source code: iscoloring.c
1: #define PETSCVEC_DLL
3: #include petscsys.h
4: #include petscis.h
6: const char *ISColoringTypes[] = {"global","ghosted","ISColoringType","IS_COLORING_",0};
10: /*@
11: ISColoringDestroy - Destroys a coloring context.
13: Collective on ISColoring
15: Input Parameter:
16: . iscoloring - the coloring context
18: Level: advanced
20: .seealso: ISColoringView(), MatGetColoring()
21: @*/
22: PetscErrorCode ISColoringDestroy(ISColoring iscoloring)
23: {
24: PetscInt i;
29: if (--iscoloring->refct > 0) return(0);
31: if (iscoloring->is) {
32: for (i=0; i<iscoloring->n; i++) {
33: ISDestroy(iscoloring->is[i]);
34: }
35: PetscFree(iscoloring->is);
36: }
37: PetscFree(iscoloring->colors);
38: PetscCommDestroy(&iscoloring->comm);
39: PetscFree(iscoloring);
40: return(0);
41: }
45: /*@C
46: ISColoringView - Views a coloring context.
48: Collective on ISColoring
50: Input Parameters:
51: + iscoloring - the coloring context
52: - viewer - the viewer
54: Level: advanced
56: .seealso: ISColoringDestroy(), ISColoringGetIS(), MatGetColoring()
57: @*/
58: PetscErrorCode ISColoringView(ISColoring iscoloring,PetscViewer viewer)
59: {
60: PetscInt i;
62: PetscTruth iascii;
63: IS *is;
67: if (!viewer) {
68: PetscViewerASCIIGetStdout(iscoloring->comm,&viewer);
69: }
72: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
73: if (iascii) {
74: MPI_Comm comm;
75: PetscMPIInt rank;
76: PetscObjectGetComm((PetscObject)viewer,&comm);
77: MPI_Comm_rank(comm,&rank);
78: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of colors %d\n",rank,iscoloring->n);
79: PetscViewerFlush(viewer);
80: } else {
81: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for ISColoring",((PetscObject)viewer)->type_name);
82: }
84: ISColoringGetIS(iscoloring,PETSC_IGNORE,&is);
85: for (i=0; i<iscoloring->n; i++) {
86: ISView(iscoloring->is[i],viewer);
87: }
88: ISColoringRestoreIS(iscoloring,&is);
89: return(0);
90: }
94: /*@C
95: ISColoringGetIS - Extracts index sets from the coloring context
97: Collective on ISColoring
99: Input Parameter:
100: . iscoloring - the coloring context
102: Output Parameters:
103: + nn - number of index sets in the coloring context
104: - is - array of index sets
106: Level: advanced
108: .seealso: ISColoringRestoreIS(), ISColoringView()
109: @*/
110: PetscErrorCode ISColoringGetIS(ISColoring iscoloring,PetscInt *nn,IS *isis[])
111: {
117: if (nn) *nn = iscoloring->n;
118: if (isis) {
119: if (!iscoloring->is) {
120: PetscInt *mcolors,**ii,nc = iscoloring->n,i,base, n = iscoloring->N;
121: ISColoringValue *colors = iscoloring->colors;
122: IS *is;
124: #if defined(PETSC_USE_DEBUG)
125: for (i=0; i<n; i++) {
126: if (((PetscInt)colors[i]) >= nc) {
127: SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"Coloring is our of range index %d value %d number colors %d",(int)i,(int)colors[i],(int)nc);
128: }
129: }
130: #endif
131:
132: /* generate the lists of nodes for each color */
133: PetscMalloc(nc*sizeof(PetscInt),&mcolors);
134: PetscMemzero(mcolors,nc*sizeof(PetscInt));
135: for (i=0; i<n; i++) {
136: mcolors[colors[i]]++;
137: }
139: PetscMalloc(nc*sizeof(PetscInt*),&ii);
140: PetscMalloc(n*sizeof(PetscInt),&ii[0]);
141: for (i=1; i<nc; i++) {
142: ii[i] = ii[i-1] + mcolors[i-1];
143: }
144: PetscMemzero(mcolors,nc*sizeof(PetscInt));
146: if (iscoloring->ctype == IS_COLORING_GLOBAL){
147: MPI_Scan(&iscoloring->N,&base,1,MPIU_INT,MPI_SUM,iscoloring->comm);
148: base -= iscoloring->N;
149: for (i=0; i<n; i++) {
150: ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
151: }
152: } else if (iscoloring->ctype == IS_COLORING_GHOSTED){
153: for (i=0; i<n; i++) {
154: ii[colors[i]][mcolors[colors[i]]++] = i; /* local idx */
155: }
156: } else {
157: SETERRQ(PETSC_ERR_SUP,"Not provided for this ISColoringType type");
158: }
159:
160: PetscMalloc(nc*sizeof(IS),&is);
161: for (i=0; i<nc; i++) {
162: ISCreateGeneral(iscoloring->comm,mcolors[i],ii[i],is+i);
163: }
165: iscoloring->is = is;
166: PetscFree(ii[0]);
167: PetscFree(ii);
168: PetscFree(mcolors);
169: }
170: *isis = iscoloring->is;
171: }
172: return(0);
173: }
177: /*@C
178: ISColoringRestoreIS - Restores the index sets extracted from the coloring context
180: Collective on ISColoring
182: Input Parameter:
183: + iscoloring - the coloring context
184: - is - array of index sets
186: Level: advanced
188: .seealso: ISColoringGetIS(), ISColoringView()
189: @*/
190: PetscErrorCode ISColoringRestoreIS(ISColoring iscoloring,IS *is[])
191: {
194:
195: /* currently nothing is done here */
197: return(0);
198: }
203: /*@C
204: ISColoringCreate - Generates an ISColoring context from lists (provided
205: by each processor) of colors for each node.
207: Collective on MPI_Comm
209: Input Parameters:
210: + comm - communicator for the processors creating the coloring
211: . ncolors - max color value
212: . n - number of nodes on this processor
213: - colors - array containing the colors for this processor, color
214: numbers begin at 0. In C/C++ this array must have been obtained with PetscMalloc()
215: and should NOT be freed (The ISColoringDestroy() will free it).
217: Output Parameter:
218: . iscoloring - the resulting coloring data structure
220: Options Database Key:
221: . -is_coloring_view - Activates ISColoringView()
223: Level: advanced
224:
225: Notes: By default sets coloring type to IS_COLORING_GLOBAL
227: .seealso: MatColoringCreate(), ISColoringView(), ISColoringDestroy(), ISColoringSetType()
229: @*/
230: PetscErrorCode ISColoringCreate(MPI_Comm comm,PetscInt ncolors,PetscInt n,const ISColoringValue colors[],ISColoring *iscoloring)
231: {
233: PetscMPIInt size,rank,tag;
234: PetscInt base,top,i;
235: PetscInt nc,ncwork;
236: PetscTruth flg;
237: MPI_Status status;
240: if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
241: if (ncolors > 65535) {
242: SETERRQ2(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);
243: } else {
244: SETERRQ2(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);
245: }
246: }
247: PetscNew(struct _n_ISColoring,iscoloring);
248: PetscCommDuplicate(comm,&(*iscoloring)->comm,&tag);
249: comm = (*iscoloring)->comm;
251: /* compute the number of the first node on my processor */
252: MPI_Comm_size(comm,&size);
254: /* should use MPI_Scan() */
255: MPI_Comm_rank(comm,&rank);
256: if (!rank) {
257: base = 0;
258: top = n;
259: } else {
260: MPI_Recv(&base,1,MPIU_INT,rank-1,tag,comm,&status);
261: top = base+n;
262: }
263: if (rank < size-1) {
264: MPI_Send(&top,1,MPIU_INT,rank+1,tag,comm);
265: }
267: /* compute the total number of colors */
268: ncwork = 0;
269: for (i=0; i<n; i++) {
270: if (ncwork < colors[i]) ncwork = colors[i];
271: }
272: ncwork++;
273: MPI_Allreduce(&ncwork,&nc,1,MPIU_INT,MPI_MAX,comm);
274: if (nc > ncolors) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Number of colors passed in %D is less then the actual number of colors in array %D",ncolors,nc);
275: (*iscoloring)->n = nc;
276: (*iscoloring)->is = 0;
277: (*iscoloring)->colors = (ISColoringValue *)colors;
278: (*iscoloring)->N = n;
279: (*iscoloring)->refct = 1;
280: (*iscoloring)->ctype = IS_COLORING_GLOBAL;
282: PetscOptionsHasName(PETSC_NULL,"-is_coloring_view",&flg);
283: if (flg) {
284: PetscViewer viewer;
285: PetscViewerASCIIGetStdout((*iscoloring)->comm,&viewer);
286: ISColoringView(*iscoloring,viewer);
287: }
288: PetscInfo1(0,"Number of colors %d\n",nc);
289: return(0);
290: }
294: /*@
295: ISPartitioningToNumbering - Takes an ISPartitioning and on each processor
296: generates an IS that contains a new global node number for each index based
297: on the partitioing.
299: Collective on IS
301: Input Parameters
302: . partitioning - a partitioning as generated by MatPartitioningApply()
304: Output Parameter:
305: . is - on each processor the index set that defines the global numbers
306: (in the new numbering) for all the nodes currently (before the partitioning)
307: on that processor
309: Level: advanced
311: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningCount()
313: @*/
314: PetscErrorCode ISPartitioningToNumbering(IS part,IS *is)
315: {
316: MPI_Comm comm;
317: PetscInt i,*indices = PETSC_NULL,np,npt,n,*starts = PETSC_NULL,*sums = PETSC_NULL,*lsizes = PETSC_NULL,*newi = PETSC_NULL;
321: PetscObjectGetComm((PetscObject)part,&comm);
323: /* count the number of partitions, i.e., virtual processors */
324: ISGetLocalSize(part,&n);
325: ISGetIndices(part,&indices);
326: np = 0;
327: for (i=0; i<n; i++) {
328: np = PetscMax(np,indices[i]);
329: }
330: MPI_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
331: np = npt+1; /* so that it looks like a MPI_Comm_size output */
333: /*
334: lsizes - number of elements of each partition on this particular processor
335: sums - total number of "previous" nodes for any particular partition
336: starts - global number of first element in each partition on this processor
337: */
338: PetscMalloc3(np,PetscInt,&lsizes,np,PetscInt,&starts,np,PetscInt,&sums);
339: PetscMemzero(lsizes,np*sizeof(PetscInt));
340: for (i=0; i<n; i++) {
341: lsizes[indices[i]]++;
342: }
343: MPI_Allreduce(lsizes,sums,np,MPIU_INT,MPI_SUM,comm);
344: MPI_Scan(lsizes,starts,np,MPIU_INT,MPI_SUM,comm);
345: for (i=0; i<np; i++) {
346: starts[i] -= lsizes[i];
347: }
348: for (i=1; i<np; i++) {
349: sums[i] += sums[i-1];
350: starts[i] += sums[i-1];
351: }
353: /*
354: For each local index give it the new global number
355: */
356: PetscMalloc(n*sizeof(PetscInt),&newi);
357: for (i=0; i<n; i++) {
358: newi[i] = starts[indices[i]]++;
359: }
360: PetscFree3(lsizes,starts,sums);
362: ISRestoreIndices(part,&indices);
363: ISCreateGeneral(comm,n,newi,is);
364: PetscFree(newi);
365: ISSetPermutation(*is);
366: return(0);
367: }
371: /*@
372: ISPartitioningCount - Takes a ISPartitioning and determines the number of
373: resulting elements on each processor
375: Collective on IS
377: Input Parameters:
378: . partitioning - a partitioning as generated by MatPartitioningApply()
380: Output Parameter:
381: . count - array of length size, to contain the number of elements assigned
382: to each partition, where size is the number of partitions generated
383: (see notes below).
385: Level: advanced
387: Notes:
388: By default the number of partitions generated (and thus the length
389: of count) is the size of the communicator associated with IS,
390: but it can be set by MatPartitioningSetNParts. The resulting array
391: of lengths can for instance serve as input of PCBJacobiSetTotalBlocks.
394: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningToNumbering(),
395: MatPartitioningSetNParts()
397: @*/
398: PetscErrorCode ISPartitioningCount(IS part,PetscInt count[])
399: {
400: MPI_Comm comm;
401: PetscInt i,*indices,np,npt,n,*lsizes;
405: PetscObjectGetComm((PetscObject)part,&comm);
407: /* count the number of partitions */
408: ISGetLocalSize(part,&n);
409: ISGetIndices(part,&indices);
410: np = 0;
411: for (i=0; i<n; i++) {
412: np = PetscMax(np,indices[i]);
413: }
414: MPI_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
415: np = npt+1; /* so that it looks like a MPI_Comm_size output */
417: /*
418: lsizes - number of elements of each partition on this particular processor
419: sums - total number of "previous" nodes for any particular partition
420: starts - global number of first element in each partition on this processor
421: */
422: PetscMalloc(np*sizeof(PetscInt),&lsizes);
423: PetscMemzero(lsizes,np*sizeof(PetscInt));
424: for (i=0; i<n; i++) {
425: lsizes[indices[i]]++;
426: }
427: ISRestoreIndices(part,&indices);
428: MPI_Allreduce(lsizes,count,np,MPIU_INT,MPI_SUM,comm);
429: PetscFree(lsizes);
430: return(0);
431: }
435: /*@
436: ISAllGather - Given an index set (IS) on each processor, generates a large
437: index set (same on each processor) by concatenating together each
438: processors index set.
440: Collective on IS
442: Input Parameter:
443: . is - the distributed index set
445: Output Parameter:
446: . isout - the concatenated index set (same on all processors)
448: Notes:
449: ISAllGather() is clearly not scalable for large index sets.
451: The IS created on each processor must be created with a common
452: communicator (e.g., PETSC_COMM_WORLD). If the index sets were created
453: with PETSC_COMM_SELF, this routine will not work as expected, since
454: each process will generate its own new IS that consists only of
455: itself.
457: Level: intermediate
459: Concepts: gather^index sets
460: Concepts: index sets^gathering to all processors
461: Concepts: IS^gathering to all processors
463: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGatherIndices()
464: @*/
465: PetscErrorCode ISAllGather(IS is,IS *isout)
466: {
468: PetscInt *indices,n,*lindices,i,N,step,first;
469: MPI_Comm comm;
470: PetscMPIInt size,*sizes = PETSC_NULL,*offsets = PETSC_NULL,nn;
471: PetscTruth stride;
477: PetscObjectGetComm((PetscObject)is,&comm);
478: MPI_Comm_size(comm,&size);
479: ISGetLocalSize(is,&n);
480: ISStride(is,&stride);
481: if (size == 1 && stride) { /* should handle parallel ISStride also */
482: ISStrideGetInfo(is,&first,&step);
483: ISCreateStride(PETSC_COMM_SELF,n,first,step,isout);
484: } else {
485: PetscMalloc2(size,PetscMPIInt,&sizes,size,PetscMPIInt,&offsets);
486:
487: nn = (PetscMPIInt)n;
488: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
489: offsets[0] = 0;
490: for (i=1;i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
491: N = offsets[size-1] + sizes[size-1];
492:
493: PetscMalloc(N*sizeof(PetscInt),&indices);
494: ISGetIndices(is,&lindices);
495: MPI_Allgatherv(lindices,nn,MPIU_INT,indices,sizes,offsets,MPIU_INT,comm);
496: ISRestoreIndices(is,&lindices);
497: PetscFree(sizes);
499: ISCreateGeneral(PETSC_COMM_SELF,N,indices,isout);
500: PetscFree2(indices,offsets);
501: }
502: return(0);
503: }
507: /*@C
508: ISAllGatherIndices - Given a a set of integers on each processor, generates a large
509: set (same on each processor) by concatenating together each processors integers
511: Collective on MPI_Comm
513: Input Parameter:
514: + comm - communicator to share the indices
515: . n - local size of set
516: - lindices - local indices
518: Output Parameter:
519: + outN - total number of indices
520: - outindices - all of the integers
522: Notes:
523: ISAllGatherIndices() is clearly not scalable for large index sets.
526: Level: intermediate
528: Concepts: gather^index sets
529: Concepts: index sets^gathering to all processors
530: Concepts: IS^gathering to all processors
532: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
533: @*/
534: PetscErrorCode ISAllGatherIndices(MPI_Comm comm,PetscInt n,const PetscInt lindices[],PetscInt *outN,PetscInt *outindices[])
535: {
537: PetscInt *indices,i,N;
538: PetscMPIInt size,*sizes = PETSC_NULL,*offsets = PETSC_NULL,nn;
541: MPI_Comm_size(comm,&size);
542: PetscMalloc2(size,PetscMPIInt,&sizes,size,PetscMPIInt,&offsets);
543:
544: nn = n;
545: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
546: offsets[0] = 0;
547: for (i=1;i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
548: N = offsets[size-1] + sizes[size-1];
550: PetscMalloc(N*sizeof(PetscInt),&indices);
551: MPI_Allgatherv((void*)lindices,nn,MPIU_INT,indices,sizes,offsets,MPIU_INT,comm);
552: PetscFree2(sizes,offsets);
554: *outindices = indices;
555: if (outN) *outN = N;
556: return(0);
557: }
563: /*@C
564: ISAllGatherColors - Given a a set of colors on each processor, generates a large
565: set (same on each processor) by concatenating together each processors colors
567: Collective on MPI_Comm
569: Input Parameter:
570: + comm - communicator to share the indices
571: . n - local size of set
572: - lindices - local colors
574: Output Parameter:
575: + outN - total number of indices
576: - outindices - all of the colors
578: Notes:
579: ISAllGatherColors() is clearly not scalable for large index sets.
582: Level: intermediate
584: Concepts: gather^index sets
585: Concepts: index sets^gathering to all processors
586: Concepts: IS^gathering to all processors
588: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather(), ISAllGatherIndices()
589: @*/
590: PetscErrorCode ISAllGatherColors(MPI_Comm comm,PetscInt n,ISColoringValue *lindices,PetscInt *outN,ISColoringValue *outindices[])
591: {
592: ISColoringValue *indices;
593: PetscErrorCode ierr;
594: PetscInt i,N;
595: PetscMPIInt size,*offsets = PETSC_NULL,*sizes = PETSC_NULL, nn = n;
598: MPI_Comm_size(comm,&size);
599: PetscMalloc2(size,PetscMPIInt,&sizes,size,PetscMPIInt,&offsets);
600:
601: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
602: offsets[0] = 0;
603: for (i=1;i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
604: N = offsets[size-1] + sizes[size-1];
605: PetscFree2(sizes,offsets);
607: PetscMalloc((N+1)*sizeof(ISColoringValue),&indices);
608: MPI_Allgatherv(lindices,(PetscMPIInt)n,MPIU_COLORING_VALUE,indices,sizes,offsets,MPIU_COLORING_VALUE,comm);
610: *outindices = indices;
611: if (outN) *outN = N;
612: return(0);
613: }