Actual source code: iscoloring.c
2: #include <petsc/private/isimpl.h>
3: #include <petscviewer.h>
4: #include <petscsf.h>
6: const char *const ISColoringTypes[] = {"global","ghosted","ISColoringType","IS_COLORING_",NULL};
8: PetscErrorCode ISColoringReference(ISColoring coloring)
9: {
10: coloring->refct++;
11: return 0;
12: }
14: /*@C
16: ISColoringSetType - indicates if the coloring is for the local representation (including ghost points) or the global representation
18: Collective on coloring
20: Input Parameters:
21: + coloring - the coloring object
22: - type - either IS_COLORING_LOCAL or IS_COLORING_GLOBAL
24: Notes:
25: With IS_COLORING_LOCAL the coloring is in the numbering of the local vector, for IS_COLORING_GLOBAL it is in the number of the global vector
27: Level: intermediate
29: .seealso: MatFDColoringCreate(), ISColoring, ISColoringCreate(), IS_COLORING_LOCAL, IS_COLORING_GLOBAL, ISColoringGetType()
31: @*/
32: PetscErrorCode ISColoringSetType(ISColoring coloring,ISColoringType type)
33: {
34: coloring->ctype = type;
35: return 0;
36: }
38: /*@C
40: ISColoringGetType - gets if the coloring is for the local representation (including ghost points) or the global representation
42: Collective on coloring
44: Input Parameter:
45: . coloring - the coloring object
47: Output Parameter:
48: . type - either IS_COLORING_LOCAL or IS_COLORING_GLOBAL
50: Level: intermediate
52: .seealso: MatFDColoringCreate(), ISColoring, ISColoringCreate(), IS_COLORING_LOCAL, IS_COLORING_GLOBAL, ISColoringSetType()
54: @*/
55: PetscErrorCode ISColoringGetType(ISColoring coloring,ISColoringType *type)
56: {
57: *type = coloring->ctype;
58: return 0;
59: }
61: /*@
62: ISColoringDestroy - Destroys a coloring context.
64: Collective on ISColoring
66: Input Parameter:
67: . iscoloring - the coloring context
69: Level: advanced
71: .seealso: ISColoringView(), MatColoring
72: @*/
73: PetscErrorCode ISColoringDestroy(ISColoring *iscoloring)
74: {
75: PetscInt i;
77: if (!*iscoloring) return 0;
79: if (--(*iscoloring)->refct > 0) {*iscoloring = NULL; return 0;}
81: if ((*iscoloring)->is) {
82: for (i=0; i<(*iscoloring)->n; i++) {
83: ISDestroy(&(*iscoloring)->is[i]);
84: }
85: PetscFree((*iscoloring)->is);
86: }
87: if ((*iscoloring)->allocated) PetscFree((*iscoloring)->colors);
88: PetscCommDestroy(&(*iscoloring)->comm);
89: PetscFree((*iscoloring));
90: return 0;
91: }
93: /*
94: ISColoringViewFromOptions - Processes command line options to determine if/how an ISColoring object is to be viewed.
96: Collective on ISColoring
98: Input Parameters:
99: + obj - the ISColoring object
100: . prefix - prefix to use for viewing, or NULL to use prefix of 'mat'
101: - optionname - option to activate viewing
103: Level: intermediate
105: Developer Note: This cannot use PetscObjectViewFromOptions() because ISColoring is not a PetscObject
107: */
108: PetscErrorCode ISColoringViewFromOptions(ISColoring obj,PetscObject bobj,const char optionname[])
109: {
110: PetscViewer viewer;
111: PetscBool flg;
112: PetscViewerFormat format;
113: char *prefix;
115: prefix = bobj ? bobj->prefix : NULL;
116: PetscOptionsGetViewer(obj->comm,NULL,prefix,optionname,&viewer,&format,&flg);
117: if (flg) {
118: PetscViewerPushFormat(viewer,format);
119: ISColoringView(obj,viewer);
120: PetscViewerPopFormat(viewer);
121: PetscViewerDestroy(&viewer);
122: }
123: return 0;
124: }
126: /*@C
127: ISColoringView - Views a coloring context.
129: Collective on ISColoring
131: Input Parameters:
132: + iscoloring - the coloring context
133: - viewer - the viewer
135: Level: advanced
137: .seealso: ISColoringDestroy(), ISColoringGetIS(), MatColoring
138: @*/
139: PetscErrorCode ISColoringView(ISColoring iscoloring,PetscViewer viewer)
140: {
141: PetscInt i;
142: PetscBool iascii;
143: IS *is;
146: if (!viewer) {
147: PetscViewerASCIIGetStdout(iscoloring->comm,&viewer);
148: }
151: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
152: if (iascii) {
153: MPI_Comm comm;
154: PetscMPIInt size,rank;
156: PetscObjectGetComm((PetscObject)viewer,&comm);
157: MPI_Comm_size(comm,&size);
158: MPI_Comm_rank(comm,&rank);
159: PetscViewerASCIIPrintf(viewer,"ISColoring Object: %d MPI processes\n",size);
160: PetscViewerASCIIPrintf(viewer,"ISColoringType: %s\n",ISColoringTypes[iscoloring->ctype]);
161: PetscViewerASCIIPushSynchronized(viewer);
162: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of colors %" PetscInt_FMT "\n",rank,iscoloring->n);
163: PetscViewerFlush(viewer);
164: PetscViewerASCIIPopSynchronized(viewer);
165: }
167: ISColoringGetIS(iscoloring,PETSC_USE_POINTER,PETSC_IGNORE,&is);
168: for (i=0; i<iscoloring->n; i++) {
169: ISView(iscoloring->is[i],viewer);
170: }
171: ISColoringRestoreIS(iscoloring,PETSC_USE_POINTER,&is);
172: return 0;
173: }
175: /*@C
176: ISColoringGetColors - Returns an array with the color for each node
178: Not Collective
180: Input Parameter:
181: . iscoloring - the coloring context
183: Output Parameters:
184: + n - number of nodes
185: . nc - number of colors
186: - colors - color for each node
188: Level: advanced
190: .seealso: ISColoringRestoreIS(), ISColoringView(), ISColoringGetIS()
191: @*/
192: PetscErrorCode ISColoringGetColors(ISColoring iscoloring,PetscInt *n,PetscInt *nc,const ISColoringValue **colors)
193: {
196: if (n) *n = iscoloring->N;
197: if (nc) *nc = iscoloring->n;
198: if (colors) *colors = iscoloring->colors;
199: return 0;
200: }
202: /*@C
203: ISColoringGetIS - Extracts index sets from the coloring context. Each is contains the nodes of one color
205: Collective on ISColoring
207: Input Parameters:
208: + iscoloring - the coloring context
209: - mode - if this value is PETSC_OWN_POINTER then the caller owns the pointer and must free the array of IS and each IS in the array
211: Output Parameters:
212: + nn - number of index sets in the coloring context
213: - is - array of index sets
215: Level: advanced
217: .seealso: ISColoringRestoreIS(), ISColoringView(), ISColoringGetColoring()
218: @*/
219: PetscErrorCode ISColoringGetIS(ISColoring iscoloring,PetscCopyMode mode, PetscInt *nn,IS *isis[])
220: {
223: if (nn) *nn = iscoloring->n;
224: if (isis) {
225: if (!iscoloring->is) {
226: PetscInt *mcolors,**ii,nc = iscoloring->n,i,base, n = iscoloring->N;
227: ISColoringValue *colors = iscoloring->colors;
228: IS *is;
230: if (PetscDefined(USE_DEBUG)) {
231: for (i=0; i<n; i++) {
233: }
234: }
236: /* generate the lists of nodes for each color */
237: PetscCalloc1(nc,&mcolors);
238: for (i=0; i<n; i++) mcolors[colors[i]]++;
240: PetscMalloc1(nc,&ii);
241: PetscMalloc1(n,&ii[0]);
242: for (i=1; i<nc; i++) ii[i] = ii[i-1] + mcolors[i-1];
243: PetscArrayzero(mcolors,nc);
245: if (iscoloring->ctype == IS_COLORING_GLOBAL) {
246: MPI_Scan(&iscoloring->N,&base,1,MPIU_INT,MPI_SUM,iscoloring->comm);
247: base -= iscoloring->N;
248: for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
249: } else if (iscoloring->ctype == IS_COLORING_LOCAL) {
250: for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i; /* local idx */
251: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this ISColoringType type");
253: PetscMalloc1(nc,&is);
254: for (i=0; i<nc; i++) {
255: ISCreateGeneral(iscoloring->comm,mcolors[i],ii[i],PETSC_COPY_VALUES,is+i);
256: }
258: if (mode != PETSC_OWN_POINTER) iscoloring->is = is;
259: *isis = is;
260: PetscFree(ii[0]);
261: PetscFree(ii);
262: PetscFree(mcolors);
263: } else {
264: *isis = iscoloring->is;
265: }
266: }
267: return 0;
268: }
270: /*@C
271: ISColoringRestoreIS - Restores the index sets extracted from the coloring context
273: Collective on ISColoring
275: Input Parameters:
276: + iscoloring - the coloring context
277: . mode - who retains ownership of the is
278: - is - array of index sets
280: Level: advanced
282: .seealso: ISColoringGetIS(), ISColoringView()
283: @*/
284: PetscErrorCode ISColoringRestoreIS(ISColoring iscoloring,PetscCopyMode mode,IS *is[])
285: {
288: /* currently nothing is done here */
289: return 0;
290: }
292: /*@
293: ISColoringCreate - Generates an ISColoring context from lists (provided
294: by each processor) of colors for each node.
296: Collective
298: Input Parameters:
299: + comm - communicator for the processors creating the coloring
300: . ncolors - max color value
301: . n - number of nodes on this processor
302: . colors - array containing the colors for this processor, color numbers begin at 0.
303: - mode - see PetscCopyMode for meaning of this flag.
305: Output Parameter:
306: . iscoloring - the resulting coloring data structure
308: Options Database Key:
309: . -is_coloring_view - Activates ISColoringView()
311: Level: advanced
313: Notes:
314: By default sets coloring type to IS_COLORING_GLOBAL
316: .seealso: MatColoringCreate(), ISColoringView(), ISColoringDestroy(), ISColoringSetType()
318: @*/
319: PetscErrorCode ISColoringCreate(MPI_Comm comm,PetscInt ncolors,PetscInt n,const ISColoringValue colors[],PetscCopyMode mode,ISColoring *iscoloring)
320: {
321: PetscMPIInt size,rank,tag;
322: PetscInt base,top,i;
323: PetscInt nc,ncwork;
324: MPI_Status status;
326: if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
328: else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Max color value exceeds limit. Perhaps reconfigure PETSc with --with-is-color-value-type=short?\n Current max: %d user requested: %" PetscInt_FMT,PETSC_IS_COLORING_MAX,ncolors);
329: }
330: PetscNew(iscoloring);
331: PetscCommDuplicate(comm,&(*iscoloring)->comm,&tag);
332: comm = (*iscoloring)->comm;
334: /* compute the number of the first node on my processor */
335: MPI_Comm_size(comm,&size);
337: /* should use MPI_Scan() */
338: MPI_Comm_rank(comm,&rank);
339: if (rank == 0) {
340: base = 0;
341: top = n;
342: } else {
343: MPI_Recv(&base,1,MPIU_INT,rank-1,tag,comm,&status);
344: top = base+n;
345: }
346: if (rank < size-1) {
347: MPI_Send(&top,1,MPIU_INT,rank+1,tag,comm);
348: }
350: /* compute the total number of colors */
351: ncwork = 0;
352: for (i=0; i<n; i++) {
353: if (ncwork < colors[i]) ncwork = colors[i];
354: }
355: ncwork++;
356: MPIU_Allreduce(&ncwork,&nc,1,MPIU_INT,MPI_MAX,comm);
358: (*iscoloring)->n = nc;
359: (*iscoloring)->is = NULL;
360: (*iscoloring)->N = n;
361: (*iscoloring)->refct = 1;
362: (*iscoloring)->ctype = IS_COLORING_GLOBAL;
363: if (mode == PETSC_COPY_VALUES) {
364: PetscMalloc1(n,&(*iscoloring)->colors);
365: PetscLogObjectMemory((PetscObject)(*iscoloring),n*sizeof(ISColoringValue));
366: PetscArraycpy((*iscoloring)->colors,colors,n);
367: (*iscoloring)->allocated = PETSC_TRUE;
368: } else if (mode == PETSC_OWN_POINTER) {
369: (*iscoloring)->colors = (ISColoringValue*)colors;
370: (*iscoloring)->allocated = PETSC_TRUE;
371: } else {
372: (*iscoloring)->colors = (ISColoringValue*)colors;
373: (*iscoloring)->allocated = PETSC_FALSE;
374: }
375: ISColoringViewFromOptions(*iscoloring,NULL,"-is_coloring_view");
376: PetscInfo(0,"Number of colors %" PetscInt_FMT "\n",nc);
377: return 0;
378: }
380: /*@
381: ISBuildTwoSided - Takes an IS that describes where we will go. Generates an IS that contains new numbers from remote or local
382: on the IS.
384: Collective on IS
386: Input Parameters:
387: + ito - an IS describes where we will go. Negative target rank will be ignored
388: - toindx - an IS describes what indices should send. NULL means sending natural numbering
390: Output Parameter:
391: . rows - contains new numbers from remote or local
393: Level: advanced
395: .seealso: MatPartitioningCreate(), ISPartitioningToNumbering(), ISPartitioningCount()
397: @*/
398: PetscErrorCode ISBuildTwoSided(IS ito,IS toindx, IS *rows)
399: {
400: const PetscInt *ito_indices,*toindx_indices;
401: PetscInt *send_indices,rstart,*recv_indices,nrecvs,nsends;
402: PetscInt *tosizes,*fromsizes,i,j,*tosizes_tmp,*tooffsets_tmp,ito_ln;
403: PetscMPIInt *toranks,*fromranks,size,target_rank,*fromperm_newtoold,nto,nfrom;
404: PetscLayout isrmap;
405: MPI_Comm comm;
406: PetscSF sf;
407: PetscSFNode *iremote;
409: PetscObjectGetComm((PetscObject)ito,&comm);
410: MPI_Comm_size(comm,&size);
411: ISGetLocalSize(ito,&ito_ln);
412: ISGetLayout(ito,&isrmap);
413: PetscLayoutGetRange(isrmap,&rstart,NULL);
414: ISGetIndices(ito,&ito_indices);
415: PetscCalloc2(size,&tosizes_tmp,size+1,&tooffsets_tmp);
416: for (i=0; i<ito_ln; i++) {
417: if (ito_indices[i]<0) continue;
419: tosizes_tmp[ito_indices[i]]++;
420: }
421: nto = 0;
422: for (i=0; i<size; i++) {
423: tooffsets_tmp[i+1] = tooffsets_tmp[i]+tosizes_tmp[i];
424: if (tosizes_tmp[i]>0) nto++;
425: }
426: PetscCalloc2(nto,&toranks,2*nto,&tosizes);
427: nto = 0;
428: for (i=0; i<size; i++) {
429: if (tosizes_tmp[i]>0) {
430: toranks[nto] = i;
431: tosizes[2*nto] = tosizes_tmp[i];/* size */
432: tosizes[2*nto+1] = tooffsets_tmp[i];/* offset */
433: nto++;
434: }
435: }
436: nsends = tooffsets_tmp[size];
437: PetscCalloc1(nsends,&send_indices);
438: if (toindx) {
439: ISGetIndices(toindx,&toindx_indices);
440: }
441: for (i=0; i<ito_ln; i++) {
442: if (ito_indices[i]<0) continue;
443: target_rank = ito_indices[i];
444: send_indices[tooffsets_tmp[target_rank]] = toindx? toindx_indices[i]:(i+rstart);
445: tooffsets_tmp[target_rank]++;
446: }
447: if (toindx) {
448: ISRestoreIndices(toindx,&toindx_indices);
449: }
450: ISRestoreIndices(ito,&ito_indices);
451: PetscFree2(tosizes_tmp,tooffsets_tmp);
452: PetscCommBuildTwoSided(comm,2,MPIU_INT,nto,toranks,tosizes,&nfrom,&fromranks,&fromsizes);
453: PetscFree2(toranks,tosizes);
454: PetscMalloc1(nfrom,&fromperm_newtoold);
455: for (i=0; i<nfrom; i++) fromperm_newtoold[i] = i;
456: PetscSortMPIIntWithArray(nfrom,fromranks,fromperm_newtoold);
457: nrecvs = 0;
458: for (i=0; i<nfrom; i++) nrecvs += fromsizes[i*2];
459: PetscCalloc1(nrecvs,&recv_indices);
460: PetscMalloc1(nrecvs,&iremote);
461: nrecvs = 0;
462: for (i=0; i<nfrom; i++) {
463: for (j=0; j<fromsizes[2*fromperm_newtoold[i]]; j++) {
464: iremote[nrecvs].rank = fromranks[i];
465: iremote[nrecvs++].index = fromsizes[2*fromperm_newtoold[i]+1]+j;
466: }
467: }
468: PetscSFCreate(comm,&sf);
469: PetscSFSetGraph(sf,nsends,nrecvs,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
470: PetscSFSetType(sf,PETSCSFBASIC);
471: /* how to put a prefix ? */
472: PetscSFSetFromOptions(sf);
473: PetscSFBcastBegin(sf,MPIU_INT,send_indices,recv_indices,MPI_REPLACE);
474: PetscSFBcastEnd(sf,MPIU_INT,send_indices,recv_indices,MPI_REPLACE);
475: PetscSFDestroy(&sf);
476: PetscFree(fromranks);
477: PetscFree(fromsizes);
478: PetscFree(fromperm_newtoold);
479: PetscFree(send_indices);
480: if (rows) {
481: PetscSortInt(nrecvs,recv_indices);
482: ISCreateGeneral(comm,nrecvs,recv_indices,PETSC_OWN_POINTER,rows);
483: } else {
484: PetscFree(recv_indices);
485: }
486: return 0;
487: }
489: /*@
490: ISPartitioningToNumbering - Takes an ISPartitioning and on each processor
491: generates an IS that contains a new global node number for each index based
492: on the partitioing.
494: Collective on IS
496: Input Parameters:
497: . partitioning - a partitioning as generated by MatPartitioningApply()
498: or MatPartitioningApplyND()
500: Output Parameter:
501: . is - on each processor the index set that defines the global numbers
502: (in the new numbering) for all the nodes currently (before the partitioning)
503: on that processor
505: Level: advanced
507: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningCount()
509: @*/
510: PetscErrorCode ISPartitioningToNumbering(IS part,IS *is)
511: {
512: MPI_Comm comm;
513: IS ndorder;
514: PetscInt i,np,npt,n,*starts = NULL,*sums = NULL,*lsizes = NULL,*newi = NULL;
515: const PetscInt *indices = NULL;
519: /* see if the partitioning comes from nested dissection */
520: PetscObjectQuery((PetscObject)part,"_petsc_matpartitioning_ndorder",(PetscObject*)&ndorder);
521: if (ndorder) {
522: PetscObjectReference((PetscObject)ndorder);
523: *is = ndorder;
524: return 0;
525: }
527: PetscObjectGetComm((PetscObject)part,&comm);
528: /* count the number of partitions, i.e., virtual processors */
529: ISGetLocalSize(part,&n);
530: ISGetIndices(part,&indices);
531: np = 0;
532: for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
533: MPIU_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
534: np = npt+1; /* so that it looks like a MPI_Comm_size output */
536: /*
537: lsizes - number of elements of each partition on this particular processor
538: sums - total number of "previous" nodes for any particular partition
539: starts - global number of first element in each partition on this processor
540: */
541: PetscMalloc3(np,&lsizes,np,&starts,np,&sums);
542: PetscArrayzero(lsizes,np);
543: for (i=0; i<n; i++) lsizes[indices[i]]++;
544: MPIU_Allreduce(lsizes,sums,np,MPIU_INT,MPI_SUM,comm);
545: MPI_Scan(lsizes,starts,np,MPIU_INT,MPI_SUM,comm);
546: for (i=0; i<np; i++) starts[i] -= lsizes[i];
547: for (i=1; i<np; i++) {
548: sums[i] += sums[i-1];
549: starts[i] += sums[i-1];
550: }
552: /*
553: For each local index give it the new global number
554: */
555: PetscMalloc1(n,&newi);
556: for (i=0; i<n; i++) newi[i] = starts[indices[i]]++;
557: PetscFree3(lsizes,starts,sums);
559: ISRestoreIndices(part,&indices);
560: ISCreateGeneral(comm,n,newi,PETSC_OWN_POINTER,is);
561: ISSetPermutation(*is);
562: return 0;
563: }
565: /*@
566: ISPartitioningCount - Takes a ISPartitioning and determines the number of
567: resulting elements on each (partition) process
569: Collective on IS
571: Input Parameters:
572: + partitioning - a partitioning as generated by MatPartitioningApply() or
573: MatPartitioningApplyND()
574: - len - length of the array count, this is the total number of partitions
576: Output Parameter:
577: . count - array of length size, to contain the number of elements assigned
578: to each partition, where size is the number of partitions generated
579: (see notes below).
581: Level: advanced
583: Notes:
584: By default the number of partitions generated (and thus the length
585: of count) is the size of the communicator associated with IS,
586: but it can be set by MatPartitioningSetNParts. The resulting array
587: of lengths can for instance serve as input of PCBJacobiSetTotalBlocks.
588: If the partitioning has been obtained by MatPartitioningApplyND(),
589: the returned count does not include the separators.
591: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningToNumbering(),
592: MatPartitioningSetNParts(), MatPartitioningApply(), MatPartitioningApplyND()
594: @*/
595: PetscErrorCode ISPartitioningCount(IS part,PetscInt len,PetscInt count[])
596: {
597: MPI_Comm comm;
598: PetscInt i,n,*lsizes;
599: const PetscInt *indices;
600: PetscMPIInt npp;
602: PetscObjectGetComm((PetscObject)part,&comm);
603: if (len == PETSC_DEFAULT) {
604: PetscMPIInt size;
605: MPI_Comm_size(comm,&size);
606: len = (PetscInt) size;
607: }
609: /* count the number of partitions */
610: ISGetLocalSize(part,&n);
611: ISGetIndices(part,&indices);
612: if (PetscDefined(USE_DEBUG)) {
613: PetscInt np = 0,npt;
614: for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
615: MPIU_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
616: np = npt+1; /* so that it looks like a MPI_Comm_size output */
618: }
620: /*
621: lsizes - number of elements of each partition on this particular processor
622: sums - total number of "previous" nodes for any particular partition
623: starts - global number of first element in each partition on this processor
624: */
625: PetscCalloc1(len,&lsizes);
626: for (i=0; i<n; i++) {
627: if (indices[i] > -1) lsizes[indices[i]]++;
628: }
629: ISRestoreIndices(part,&indices);
630: PetscMPIIntCast(len,&npp);
631: MPIU_Allreduce(lsizes,count,npp,MPIU_INT,MPI_SUM,comm);
632: PetscFree(lsizes);
633: return 0;
634: }
636: /*@
637: ISAllGather - Given an index set (IS) on each processor, generates a large
638: index set (same on each processor) by concatenating together each
639: processors index set.
641: Collective on IS
643: Input Parameter:
644: . is - the distributed index set
646: Output Parameter:
647: . isout - the concatenated index set (same on all processors)
649: Notes:
650: ISAllGather() is clearly not scalable for large index sets.
652: The IS created on each processor must be created with a common
653: communicator (e.g., PETSC_COMM_WORLD). If the index sets were created
654: with PETSC_COMM_SELF, this routine will not work as expected, since
655: each process will generate its own new IS that consists only of
656: itself.
658: The communicator for this new IS is PETSC_COMM_SELF
660: Level: intermediate
662: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock()
663: @*/
664: PetscErrorCode ISAllGather(IS is,IS *isout)
665: {
666: PetscInt *indices,n,i,N,step,first;
667: const PetscInt *lindices;
668: MPI_Comm comm;
669: PetscMPIInt size,*sizes = NULL,*offsets = NULL,nn;
670: PetscBool stride;
675: PetscObjectGetComm((PetscObject)is,&comm);
676: MPI_Comm_size(comm,&size);
677: ISGetLocalSize(is,&n);
678: PetscObjectTypeCompare((PetscObject)is,ISSTRIDE,&stride);
679: if (size == 1 && stride) { /* should handle parallel ISStride also */
680: ISStrideGetInfo(is,&first,&step);
681: ISCreateStride(PETSC_COMM_SELF,n,first,step,isout);
682: } else {
683: PetscMalloc2(size,&sizes,size,&offsets);
685: PetscMPIIntCast(n,&nn);
686: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
687: offsets[0] = 0;
688: for (i=1; i<size; i++) {
689: PetscInt s = offsets[i-1] + sizes[i-1];
690: PetscMPIIntCast(s,&offsets[i]);
691: }
692: N = offsets[size-1] + sizes[size-1];
694: PetscMalloc1(N,&indices);
695: ISGetIndices(is,&lindices);
696: MPI_Allgatherv((void*)lindices,nn,MPIU_INT,indices,sizes,offsets,MPIU_INT,comm);
697: ISRestoreIndices(is,&lindices);
698: PetscFree2(sizes,offsets);
700: ISCreateGeneral(PETSC_COMM_SELF,N,indices,PETSC_OWN_POINTER,isout);
701: }
702: return 0;
703: }
705: /*@C
706: ISAllGatherColors - Given a a set of colors on each processor, generates a large
707: set (same on each processor) by concatenating together each processors colors
709: Collective
711: Input Parameters:
712: + comm - communicator to share the indices
713: . n - local size of set
714: - lindices - local colors
716: Output Parameters:
717: + outN - total number of indices
718: - outindices - all of the colors
720: Notes:
721: ISAllGatherColors() is clearly not scalable for large index sets.
723: Level: intermediate
725: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
726: @*/
727: PetscErrorCode ISAllGatherColors(MPI_Comm comm,PetscInt n,ISColoringValue *lindices,PetscInt *outN,ISColoringValue *outindices[])
728: {
729: ISColoringValue *indices;
730: PetscInt i,N;
731: PetscMPIInt size,*offsets = NULL,*sizes = NULL, nn = n;
733: MPI_Comm_size(comm,&size);
734: PetscMalloc2(size,&sizes,size,&offsets);
736: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
737: offsets[0] = 0;
738: for (i=1; i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
739: N = offsets[size-1] + sizes[size-1];
740: PetscFree2(sizes,offsets);
742: PetscMalloc1(N+1,&indices);
743: MPI_Allgatherv(lindices,(PetscMPIInt)n,MPIU_COLORING_VALUE,indices,sizes,offsets,MPIU_COLORING_VALUE,comm);
745: *outindices = indices;
746: if (outN) *outN = N;
747: return 0;
748: }
750: /*@
751: ISComplement - Given an index set (IS) generates the complement index set. That is all
752: all indices that are NOT in the given set.
754: Collective on IS
756: Input Parameters:
757: + is - the index set
758: . nmin - the first index desired in the local part of the complement
759: - 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)
761: Output Parameter:
762: . isout - the complement
764: Notes:
765: The communicator for this new IS is the same as for the input IS
767: For a parallel IS, this will generate the local part of the complement on each process
769: To generate the entire complement (on each process) of a parallel IS, first call ISAllGather() and then
770: call this routine.
772: Level: intermediate
774: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
775: @*/
776: PetscErrorCode ISComplement(IS is,PetscInt nmin,PetscInt nmax,IS *isout)
777: {
778: const PetscInt *indices;
779: PetscInt n,i,j,unique,cnt,*nindices;
780: PetscBool sorted;
786: ISSorted(is,&sorted);
789: ISGetLocalSize(is,&n);
790: ISGetIndices(is,&indices);
791: if (PetscDefined(USE_DEBUG)) {
792: for (i=0; i<n; i++) {
795: }
796: }
797: /* Count number of unique entries */
798: unique = (n>0);
799: for (i=0; i<n-1; i++) {
800: if (indices[i+1] != indices[i]) unique++;
801: }
802: PetscMalloc1(nmax-nmin-unique,&nindices);
803: cnt = 0;
804: for (i=nmin,j=0; i<nmax; i++) {
805: if (j<n && i==indices[j]) do { j++; } while (j<n && i==indices[j]);
806: else nindices[cnt++] = i;
807: }
809: ISCreateGeneral(PetscObjectComm((PetscObject)is),cnt,nindices,PETSC_OWN_POINTER,isout);
810: ISRestoreIndices(is,&indices);
811: return 0;
812: }