Actual source code: iscoloring.c
petsc-3.6.4 2016-04-12
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(), MatColoring
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: if ((*iscoloring)->allocated) {PetscFree((*iscoloring)->colors);}
56: PetscCommDestroy(&(*iscoloring)->comm);
57: PetscFree((*iscoloring));
58: return(0);
59: }
63: /*
64: ISColoringViewFromOptions - Processes command line options to determine if/how an ISColoring object is to be viewed.
66: Collective on ISColoring
68: Input Parameters:
69: + obj - the ISColoring object
70: . prefix - prefix to use for viewing, or NULL to use prefix of 'mat'
71: - optionname - option to activate viewing
73: Level: intermediate
75: Developer Note: This cannot use PetscObjectViewFromOptions() because ISColoring is not a PetscObject
77: */
78: PetscErrorCode ISColoringViewFromOptions(ISColoring obj,PetscObject bobj,const char optionname[])
79: {
80: PetscErrorCode ierr;
81: PetscViewer viewer;
82: PetscBool flg;
83: PetscViewerFormat format;
84: char *prefix;
87: prefix = bobj ? bobj->prefix : NULL;
88: PetscOptionsGetViewer(obj->comm,prefix,optionname,&viewer,&format,&flg);
89: if (flg) {
90: PetscViewerPushFormat(viewer,format);
91: ISColoringView(obj,viewer);
92: PetscViewerPopFormat(viewer);
93: PetscViewerDestroy(&viewer);
94: }
95: return(0);
96: }
100: /*@C
101: ISColoringView - Views a coloring context.
103: Collective on ISColoring
105: Input Parameters:
106: + iscoloring - the coloring context
107: - viewer - the viewer
109: Level: advanced
111: .seealso: ISColoringDestroy(), ISColoringGetIS(), MatColoring
112: @*/
113: PetscErrorCode ISColoringView(ISColoring iscoloring,PetscViewer viewer)
114: {
115: PetscInt i;
117: PetscBool iascii;
118: IS *is;
122: if (!viewer) {
123: PetscViewerASCIIGetStdout(iscoloring->comm,&viewer);
124: }
127: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
128: if (iascii) {
129: MPI_Comm comm;
130: PetscMPIInt size,rank;
132: PetscObjectGetComm((PetscObject)viewer,&comm);
133: MPI_Comm_size(comm,&size);
134: MPI_Comm_rank(comm,&rank);
135: PetscViewerASCIIPrintf(viewer,"ISColoring Object: %d MPI processes\n",size);
136: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
137: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of colors %d\n",rank,iscoloring->n);
138: PetscViewerFlush(viewer);
139: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
140: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for ISColoring",((PetscObject)viewer)->type_name);
142: ISColoringGetIS(iscoloring,PETSC_IGNORE,&is);
143: for (i=0; i<iscoloring->n; i++) {
144: ISView(iscoloring->is[i],viewer);
145: }
146: ISColoringRestoreIS(iscoloring,&is);
147: return(0);
148: }
152: /*@C
153: ISColoringGetIS - Extracts index sets from the coloring context
155: Collective on ISColoring
157: Input Parameter:
158: . iscoloring - the coloring context
160: Output Parameters:
161: + nn - number of index sets in the coloring context
162: - is - array of index sets
164: Level: advanced
166: .seealso: ISColoringRestoreIS(), ISColoringView()
167: @*/
168: PetscErrorCode ISColoringGetIS(ISColoring iscoloring,PetscInt *nn,IS *isis[])
169: {
175: if (nn) *nn = iscoloring->n;
176: if (isis) {
177: if (!iscoloring->is) {
178: PetscInt *mcolors,**ii,nc = iscoloring->n,i,base, n = iscoloring->N;
179: ISColoringValue *colors = iscoloring->colors;
180: IS *is;
182: #if defined(PETSC_USE_DEBUG)
183: for (i=0; i<n; i++) {
184: 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);
185: }
186: #endif
188: /* generate the lists of nodes for each color */
189: PetscCalloc1(nc,&mcolors);
190: for (i=0; i<n; i++) mcolors[colors[i]]++;
192: PetscMalloc1(nc,&ii);
193: PetscMalloc1(n,&ii[0]);
194: for (i=1; i<nc; i++) ii[i] = ii[i-1] + mcolors[i-1];
195: PetscMemzero(mcolors,nc*sizeof(PetscInt));
197: if (iscoloring->ctype == IS_COLORING_GLOBAL) {
198: MPI_Scan(&iscoloring->N,&base,1,MPIU_INT,MPI_SUM,iscoloring->comm);
199: base -= iscoloring->N;
200: for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
201: } else if (iscoloring->ctype == IS_COLORING_GHOSTED) {
202: for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i; /* local idx */
203: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this ISColoringType type");
205: PetscMalloc1(nc,&is);
206: for (i=0; i<nc; i++) {
207: ISCreateGeneral(iscoloring->comm,mcolors[i],ii[i],PETSC_COPY_VALUES,is+i);
208: }
210: iscoloring->is = is;
211: PetscFree(ii[0]);
212: PetscFree(ii);
213: PetscFree(mcolors);
214: }
215: *isis = iscoloring->is;
216: }
217: return(0);
218: }
222: /*@C
223: ISColoringRestoreIS - Restores the index sets extracted from the coloring context
225: Collective on ISColoring
227: Input Parameter:
228: + iscoloring - the coloring context
229: - is - array of index sets
231: Level: advanced
233: .seealso: ISColoringGetIS(), ISColoringView()
234: @*/
235: PetscErrorCode ISColoringRestoreIS(ISColoring iscoloring,IS *is[])
236: {
240: /* currently nothing is done here */
241: return(0);
242: }
247: /*@
248: ISColoringCreate - Generates an ISColoring context from lists (provided
249: by each processor) of colors for each node.
251: Collective on MPI_Comm
253: Input Parameters:
254: + comm - communicator for the processors creating the coloring
255: . ncolors - max color value
256: . n - number of nodes on this processor
257: . colors - array containing the colors for this processor, color numbers begin at 0.
258: - mode - see PetscCopyMode for meaning of this flag.
260: Output Parameter:
261: . iscoloring - the resulting coloring data structure
263: Options Database Key:
264: . -is_coloring_view - Activates ISColoringView()
266: Level: advanced
268: Notes: By default sets coloring type to IS_COLORING_GLOBAL
270: .seealso: MatColoringCreate(), ISColoringView(), ISColoringDestroy(), ISColoringSetType()
272: @*/
273: PetscErrorCode ISColoringCreate(MPI_Comm comm,PetscInt ncolors,PetscInt n,const ISColoringValue colors[],PetscCopyMode mode,ISColoring *iscoloring)
274: {
276: PetscMPIInt size,rank,tag;
277: PetscInt base,top,i;
278: PetscInt nc,ncwork;
279: MPI_Status status;
282: if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
283: 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);
284: 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);
285: }
286: PetscNew(iscoloring);
287: PetscCommDuplicate(comm,&(*iscoloring)->comm,&tag);
288: comm = (*iscoloring)->comm;
290: /* compute the number of the first node on my processor */
291: MPI_Comm_size(comm,&size);
293: /* should use MPI_Scan() */
294: MPI_Comm_rank(comm,&rank);
295: if (!rank) {
296: base = 0;
297: top = n;
298: } else {
299: MPI_Recv(&base,1,MPIU_INT,rank-1,tag,comm,&status);
300: top = base+n;
301: }
302: if (rank < size-1) {
303: MPI_Send(&top,1,MPIU_INT,rank+1,tag,comm);
304: }
306: /* compute the total number of colors */
307: ncwork = 0;
308: for (i=0; i<n; i++) {
309: if (ncwork < colors[i]) ncwork = colors[i];
310: }
311: ncwork++;
312: MPI_Allreduce(&ncwork,&nc,1,MPIU_INT,MPI_MAX,comm);
313: 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);
314: (*iscoloring)->n = nc;
315: (*iscoloring)->is = 0;
316: (*iscoloring)->N = n;
317: (*iscoloring)->refct = 1;
318: (*iscoloring)->ctype = IS_COLORING_GLOBAL;
319: if (mode == PETSC_COPY_VALUES) {
320: PetscMalloc1(n,&(*iscoloring)->colors);
321: PetscLogObjectMemory((PetscObject)(*iscoloring),n*sizeof(ISColoringValue));
322: PetscMemcpy((*iscoloring)->colors,colors,n*sizeof(ISColoringValue));
323: (*iscoloring)->allocated = PETSC_TRUE;
324: } else if (mode == PETSC_OWN_POINTER) {
325: (*iscoloring)->colors = (ISColoringValue*)colors;
326: (*iscoloring)->allocated = PETSC_TRUE;
327: } else {
328: (*iscoloring)->colors = (ISColoringValue*)colors;
329: (*iscoloring)->allocated = PETSC_FALSE;
330: }
331: ISColoringViewFromOptions(*iscoloring,NULL,"-is_coloring_view");
332: PetscInfo1(0,"Number of colors %D\n",nc);
333: return(0);
334: }
338: /*@
339: ISPartitioningToNumbering - Takes an ISPartitioning and on each processor
340: generates an IS that contains a new global node number for each index based
341: on the partitioing.
343: Collective on IS
345: Input Parameters
346: . partitioning - a partitioning as generated by MatPartitioningApply()
348: Output Parameter:
349: . is - on each processor the index set that defines the global numbers
350: (in the new numbering) for all the nodes currently (before the partitioning)
351: on that processor
353: Level: advanced
355: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningCount()
357: @*/
358: PetscErrorCode ISPartitioningToNumbering(IS part,IS *is)
359: {
360: MPI_Comm comm;
361: PetscInt i,np,npt,n,*starts = NULL,*sums = NULL,*lsizes = NULL,*newi = NULL;
362: const PetscInt *indices = NULL;
366: PetscObjectGetComm((PetscObject)part,&comm);
368: /* count the number of partitions, i.e., virtual processors */
369: ISGetLocalSize(part,&n);
370: ISGetIndices(part,&indices);
371: np = 0;
372: for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
373: MPI_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
374: np = npt+1; /* so that it looks like a MPI_Comm_size output */
376: /*
377: lsizes - number of elements of each partition on this particular processor
378: sums - total number of "previous" nodes for any particular partition
379: starts - global number of first element in each partition on this processor
380: */
381: PetscMalloc3(np,&lsizes,np,&starts,np,&sums);
382: PetscMemzero(lsizes,np*sizeof(PetscInt));
383: for (i=0; i<n; i++) lsizes[indices[i]]++;
384: MPI_Allreduce(lsizes,sums,np,MPIU_INT,MPI_SUM,comm);
385: MPI_Scan(lsizes,starts,np,MPIU_INT,MPI_SUM,comm);
386: for (i=0; i<np; i++) starts[i] -= lsizes[i];
387: for (i=1; i<np; i++) {
388: sums[i] += sums[i-1];
389: starts[i] += sums[i-1];
390: }
392: /*
393: For each local index give it the new global number
394: */
395: PetscMalloc1(n,&newi);
396: for (i=0; i<n; i++) newi[i] = starts[indices[i]]++;
397: PetscFree3(lsizes,starts,sums);
399: ISRestoreIndices(part,&indices);
400: ISCreateGeneral(comm,n,newi,PETSC_OWN_POINTER,is);
401: ISSetPermutation(*is);
402: return(0);
403: }
407: /*@
408: ISPartitioningCount - Takes a ISPartitioning and determines the number of
409: resulting elements on each (partition) process
411: Collective on IS
413: Input Parameters:
414: + partitioning - a partitioning as generated by MatPartitioningApply()
415: - len - length of the array count, this is the total number of partitions
417: Output Parameter:
418: . count - array of length size, to contain the number of elements assigned
419: to each partition, where size is the number of partitions generated
420: (see notes below).
422: Level: advanced
424: Notes:
425: By default the number of partitions generated (and thus the length
426: of count) is the size of the communicator associated with IS,
427: but it can be set by MatPartitioningSetNParts. The resulting array
428: of lengths can for instance serve as input of PCBJacobiSetTotalBlocks.
431: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningToNumbering(),
432: MatPartitioningSetNParts(), MatPartitioningApply()
434: @*/
435: PetscErrorCode ISPartitioningCount(IS part,PetscInt len,PetscInt count[])
436: {
437: MPI_Comm comm;
438: PetscInt i,n,*lsizes;
439: const PetscInt *indices;
441: PetscMPIInt npp;
444: PetscObjectGetComm((PetscObject)part,&comm);
445: if (len == PETSC_DEFAULT) {
446: PetscMPIInt size;
447: MPI_Comm_size(comm,&size);
448: len = (PetscInt) size;
449: }
451: /* count the number of partitions */
452: ISGetLocalSize(part,&n);
453: ISGetIndices(part,&indices);
454: #if defined(PETSC_USE_DEBUG)
455: {
456: PetscInt np = 0,npt;
457: for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
458: MPI_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
459: np = npt+1; /* so that it looks like a MPI_Comm_size output */
460: 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);
461: }
462: #endif
464: /*
465: lsizes - number of elements of each partition on this particular processor
466: sums - total number of "previous" nodes for any particular partition
467: starts - global number of first element in each partition on this processor
468: */
469: PetscCalloc1(len,&lsizes);
470: for (i=0; i<n; i++) lsizes[indices[i]]++;
471: ISRestoreIndices(part,&indices);
472: PetscMPIIntCast(len,&npp);
473: MPI_Allreduce(lsizes,count,npp,MPIU_INT,MPI_SUM,comm);
474: PetscFree(lsizes);
475: return(0);
476: }
480: /*@
481: ISAllGather - Given an index set (IS) on each processor, generates a large
482: index set (same on each processor) by concatenating together each
483: processors index set.
485: Collective on IS
487: Input Parameter:
488: . is - the distributed index set
490: Output Parameter:
491: . isout - the concatenated index set (same on all processors)
493: Notes:
494: ISAllGather() is clearly not scalable for large index sets.
496: The IS created on each processor must be created with a common
497: communicator (e.g., PETSC_COMM_WORLD). If the index sets were created
498: with PETSC_COMM_SELF, this routine will not work as expected, since
499: each process will generate its own new IS that consists only of
500: itself.
502: The communicator for this new IS is PETSC_COMM_SELF
504: Level: intermediate
506: Concepts: gather^index sets
507: Concepts: index sets^gathering to all processors
508: Concepts: IS^gathering to all processors
510: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock()
511: @*/
512: PetscErrorCode ISAllGather(IS is,IS *isout)
513: {
515: PetscInt *indices,n,i,N,step,first;
516: const PetscInt *lindices;
517: MPI_Comm comm;
518: PetscMPIInt size,*sizes = NULL,*offsets = NULL,nn;
519: PetscBool stride;
525: PetscObjectGetComm((PetscObject)is,&comm);
526: MPI_Comm_size(comm,&size);
527: ISGetLocalSize(is,&n);
528: PetscObjectTypeCompare((PetscObject)is,ISSTRIDE,&stride);
529: if (size == 1 && stride) { /* should handle parallel ISStride also */
530: ISStrideGetInfo(is,&first,&step);
531: ISCreateStride(PETSC_COMM_SELF,n,first,step,isout);
532: } else {
533: PetscMalloc2(size,&sizes,size,&offsets);
535: PetscMPIIntCast(n,&nn);
536: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
537: offsets[0] = 0;
538: for (i=1; i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
539: N = offsets[size-1] + sizes[size-1];
541: PetscMalloc1(N,&indices);
542: ISGetIndices(is,&lindices);
543: MPI_Allgatherv((void*)lindices,nn,MPIU_INT,indices,sizes,offsets,MPIU_INT,comm);
544: ISRestoreIndices(is,&lindices);
545: PetscFree2(sizes,offsets);
547: ISCreateGeneral(PETSC_COMM_SELF,N,indices,PETSC_OWN_POINTER,isout);
548: }
549: return(0);
550: }
554: /*@C
555: ISAllGatherColors - Given a a set of colors on each processor, generates a large
556: set (same on each processor) by concatenating together each processors colors
558: Collective on MPI_Comm
560: Input Parameter:
561: + comm - communicator to share the indices
562: . n - local size of set
563: - lindices - local colors
565: Output Parameter:
566: + outN - total number of indices
567: - outindices - all of the colors
569: Notes:
570: ISAllGatherColors() is clearly not scalable for large index sets.
573: Level: intermediate
575: Concepts: gather^index sets
576: Concepts: index sets^gathering to all processors
577: Concepts: IS^gathering to all processors
579: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
580: @*/
581: PetscErrorCode ISAllGatherColors(MPI_Comm comm,PetscInt n,ISColoringValue *lindices,PetscInt *outN,ISColoringValue *outindices[])
582: {
583: ISColoringValue *indices;
584: PetscErrorCode ierr;
585: PetscInt i,N;
586: PetscMPIInt size,*offsets = NULL,*sizes = NULL, nn = n;
589: MPI_Comm_size(comm,&size);
590: PetscMalloc2(size,&sizes,size,&offsets);
592: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
593: offsets[0] = 0;
594: for (i=1; i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
595: N = offsets[size-1] + sizes[size-1];
596: PetscFree2(sizes,offsets);
598: PetscMalloc1(N+1,&indices);
599: MPI_Allgatherv(lindices,(PetscMPIInt)n,MPIU_COLORING_VALUE,indices,sizes,offsets,MPIU_COLORING_VALUE,comm);
601: *outindices = indices;
602: if (outN) *outN = N;
603: return(0);
604: }
608: /*@
609: ISComplement - Given an index set (IS) generates the complement index set. That is all
610: all indices that are NOT in the given set.
612: Collective on IS
614: Input Parameter:
615: + is - the index set
616: . nmin - the first index desired in the local part of the complement
617: - 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)
619: Output Parameter:
620: . isout - the complement
622: Notes: The communicator for this new IS is the same as for the input IS
624: For a parallel IS, this will generate the local part of the complement on each process
626: To generate the entire complement (on each process) of a parallel IS, first call ISAllGather() and then
627: call this routine.
629: Level: intermediate
631: Concepts: gather^index sets
632: Concepts: index sets^gathering to all processors
633: Concepts: IS^gathering to all processors
635: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
636: @*/
637: PetscErrorCode ISComplement(IS is,PetscInt nmin,PetscInt nmax,IS *isout)
638: {
640: const PetscInt *indices;
641: PetscInt n,i,j,unique,cnt,*nindices;
642: PetscBool sorted;
647: if (nmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nmin %D cannot be negative",nmin);
648: if (nmin > nmax) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nmin %D cannot be greater than nmax %D",nmin,nmax);
649: ISSorted(is,&sorted);
650: if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Index set must be sorted");
652: ISGetLocalSize(is,&n);
653: ISGetIndices(is,&indices);
654: #if defined(PETSC_USE_DEBUG)
655: for (i=0; i<n; i++) {
656: 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);
657: 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);
658: }
659: #endif
660: /* Count number of unique entries */
661: unique = (n>0);
662: for (i=0; i<n-1; i++) {
663: if (indices[i+1] != indices[i]) unique++;
664: }
665: PetscMalloc1(nmax-nmin-unique,&nindices);
666: cnt = 0;
667: for (i=nmin,j=0; i<nmax; i++) {
668: if (j<n && i==indices[j]) do { j++; } while (j<n && i==indices[j]);
669: else nindices[cnt++] = i;
670: }
671: 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);
672: ISCreateGeneral(PetscObjectComm((PetscObject)is),cnt,nindices,PETSC_OWN_POINTER,isout);
673: ISRestoreIndices(is,&indices);
674: return(0);
675: }