Actual source code: iscoloring.c
petsc-3.12.5 2020-03-29
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: /*@C
17: ISColoringSetType - indicates if the coloring is for the local representation (including ghost points) or the global representation
19: Collective on coloring
21: Input Parameters:
22: + coloring - the coloring object
23: - type - either IS_COLORING_LOCAL or IS_COLORING_GLOBAL
25: Notes:
26: 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
28: Level: intermediate
30: .seealso: MatFDColoringCreate(), ISColoring, ISColoringCreate(), IS_COLORING_LOCAL, IS_COLORING_GLOBAL, ISColoringGetType()
32: @*/
33: PetscErrorCode ISColoringSetType(ISColoring coloring,ISColoringType type)
34: {
36: coloring->ctype = type;
37: return(0);
38: }
40: /*@C
42: ISColoringGetType - gets if the coloring is for the local representation (including ghost points) or the global representation
44: Collective on coloring
46: Input Parameter:
47: . coloring - the coloring object
49: Output Parameter:
50: . type - either IS_COLORING_LOCAL or IS_COLORING_GLOBAL
52: Level: intermediate
54: .seealso: MatFDColoringCreate(), ISColoring, ISColoringCreate(), IS_COLORING_LOCAL, IS_COLORING_GLOBAL, ISColoringSetType()
56: @*/
57: PetscErrorCode ISColoringGetType(ISColoring coloring,ISColoringType *type)
58: {
60: *type = coloring->ctype;
61: return(0);
62: }
64: /*@
65: ISColoringDestroy - Destroys a coloring context.
67: Collective on ISColoring
69: Input Parameter:
70: . iscoloring - the coloring context
72: Level: advanced
74: .seealso: ISColoringView(), MatColoring
75: @*/
76: PetscErrorCode ISColoringDestroy(ISColoring *iscoloring)
77: {
78: PetscInt i;
82: if (!*iscoloring) return(0);
84: if (--(*iscoloring)->refct > 0) {*iscoloring = 0; return(0);}
86: if ((*iscoloring)->is) {
87: for (i=0; i<(*iscoloring)->n; i++) {
88: ISDestroy(&(*iscoloring)->is[i]);
89: }
90: PetscFree((*iscoloring)->is);
91: }
92: if ((*iscoloring)->allocated) {PetscFree((*iscoloring)->colors);}
93: PetscCommDestroy(&(*iscoloring)->comm);
94: PetscFree((*iscoloring));
95: return(0);
96: }
98: /*
99: ISColoringViewFromOptions - Processes command line options to determine if/how an ISColoring object is to be viewed.
101: Collective on ISColoring
103: Input Parameters:
104: + obj - the ISColoring object
105: . prefix - prefix to use for viewing, or NULL to use prefix of 'mat'
106: - optionname - option to activate viewing
108: Level: intermediate
110: Developer Note: This cannot use PetscObjectViewFromOptions() because ISColoring is not a PetscObject
112: */
113: PetscErrorCode ISColoringViewFromOptions(ISColoring obj,PetscObject bobj,const char optionname[])
114: {
115: PetscErrorCode ierr;
116: PetscViewer viewer;
117: PetscBool flg;
118: PetscViewerFormat format;
119: char *prefix;
122: prefix = bobj ? bobj->prefix : NULL;
123: PetscOptionsGetViewer(obj->comm,NULL,prefix,optionname,&viewer,&format,&flg);
124: if (flg) {
125: PetscViewerPushFormat(viewer,format);
126: ISColoringView(obj,viewer);
127: PetscViewerPopFormat(viewer);
128: PetscViewerDestroy(&viewer);
129: }
130: return(0);
131: }
133: /*@C
134: ISColoringView - Views a coloring context.
136: Collective on ISColoring
138: Input Parameters:
139: + iscoloring - the coloring context
140: - viewer - the viewer
142: Level: advanced
144: .seealso: ISColoringDestroy(), ISColoringGetIS(), MatColoring
145: @*/
146: PetscErrorCode ISColoringView(ISColoring iscoloring,PetscViewer viewer)
147: {
148: PetscInt i;
150: PetscBool iascii;
151: IS *is;
155: if (!viewer) {
156: PetscViewerASCIIGetStdout(iscoloring->comm,&viewer);
157: }
160: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
161: if (iascii) {
162: MPI_Comm comm;
163: PetscMPIInt size,rank;
165: PetscObjectGetComm((PetscObject)viewer,&comm);
166: MPI_Comm_size(comm,&size);
167: MPI_Comm_rank(comm,&rank);
168: PetscViewerASCIIPrintf(viewer,"ISColoring Object: %d MPI processes\n",size);
169: PetscViewerASCIIPrintf(viewer,"ISColoringType: %s\n",ISColoringTypes[iscoloring->ctype]);
170: PetscViewerASCIIPushSynchronized(viewer);
171: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of colors %d\n",rank,iscoloring->n);
172: PetscViewerFlush(viewer);
173: PetscViewerASCIIPopSynchronized(viewer);
174: }
176: ISColoringGetIS(iscoloring,PETSC_USE_POINTER,PETSC_IGNORE,&is);
177: for (i=0; i<iscoloring->n; i++) {
178: ISView(iscoloring->is[i],viewer);
179: }
180: ISColoringRestoreIS(iscoloring,PETSC_USE_POINTER,&is);
181: return(0);
182: }
184: /*@C
185: ISColoringGetColors - Returns an array with the color for each node
187: Not Collective
189: Input Parameter:
190: . iscoloring - the coloring context
192: Output Parameters:
193: + n - number of nodes
194: . nc - number of colors
195: - colors - color for each node
197: Level: advanced
199: .seealso: ISColoringRestoreIS(), ISColoringView(), ISColoringGetIS()
200: @*/
201: PetscErrorCode ISColoringGetColors(ISColoring iscoloring,PetscInt *n,PetscInt *nc,const ISColoringValue **colors)
202: {
206: if (n) *n = iscoloring->N;
207: if (nc) *nc = iscoloring->n;
208: if (colors) *colors = iscoloring->colors;
209: return(0);
210: }
212: /*@C
213: ISColoringGetIS - Extracts index sets from the coloring context. Each is contains the nodes of one color
215: Collective on ISColoring
217: Input Parameter:
218: + iscoloring - the coloring context
219: - 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
221: Output Parameters:
222: + nn - number of index sets in the coloring context
223: - is - array of index sets
225: Level: advanced
227: .seealso: ISColoringRestoreIS(), ISColoringView(), ISColoringGetColoring()
228: @*/
229: PetscErrorCode ISColoringGetIS(ISColoring iscoloring,PetscCopyMode mode, PetscInt *nn,IS *isis[])
230: {
236: if (nn) *nn = iscoloring->n;
237: if (isis) {
238: if (!iscoloring->is) {
239: PetscInt *mcolors,**ii,nc = iscoloring->n,i,base, n = iscoloring->N;
240: ISColoringValue *colors = iscoloring->colors;
241: IS *is;
243: #if defined(PETSC_USE_DEBUG)
244: for (i=0; i<n; i++) {
245: 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);
246: }
247: #endif
249: /* generate the lists of nodes for each color */
250: PetscCalloc1(nc,&mcolors);
251: for (i=0; i<n; i++) mcolors[colors[i]]++;
253: PetscMalloc1(nc,&ii);
254: PetscMalloc1(n,&ii[0]);
255: for (i=1; i<nc; i++) ii[i] = ii[i-1] + mcolors[i-1];
256: PetscArrayzero(mcolors,nc);
258: if (iscoloring->ctype == IS_COLORING_GLOBAL) {
259: MPI_Scan(&iscoloring->N,&base,1,MPIU_INT,MPI_SUM,iscoloring->comm);
260: base -= iscoloring->N;
261: for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
262: } else if (iscoloring->ctype == IS_COLORING_LOCAL) {
263: for (i=0; i<n; i++) ii[colors[i]][mcolors[colors[i]]++] = i; /* local idx */
264: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this ISColoringType type");
266: PetscMalloc1(nc,&is);
267: for (i=0; i<nc; i++) {
268: ISCreateGeneral(iscoloring->comm,mcolors[i],ii[i],PETSC_COPY_VALUES,is+i);
269: }
271: if (mode != PETSC_OWN_POINTER) iscoloring->is = is;
272: *isis = is;
273: PetscFree(ii[0]);
274: PetscFree(ii);
275: PetscFree(mcolors);
276: } else {
277: *isis = iscoloring->is;
278: }
279: }
280: return(0);
281: }
283: /*@C
284: ISColoringRestoreIS - Restores the index sets extracted from the coloring context
286: Collective on ISColoring
288: Input Parameter:
289: + iscoloring - the coloring context
290: . mode - who retains ownership of the is
291: - is - array of index sets
293: Level: advanced
295: .seealso: ISColoringGetIS(), ISColoringView()
296: @*/
297: PetscErrorCode ISColoringRestoreIS(ISColoring iscoloring,PetscCopyMode mode,IS *is[])
298: {
302: /* currently nothing is done here */
303: return(0);
304: }
307: /*@
308: ISColoringCreate - Generates an ISColoring context from lists (provided
309: by each processor) of colors for each node.
311: Collective
313: Input Parameters:
314: + comm - communicator for the processors creating the coloring
315: . ncolors - max color value
316: . n - number of nodes on this processor
317: . colors - array containing the colors for this processor, color numbers begin at 0.
318: - mode - see PetscCopyMode for meaning of this flag.
320: Output Parameter:
321: . iscoloring - the resulting coloring data structure
323: Options Database Key:
324: . -is_coloring_view - Activates ISColoringView()
326: Level: advanced
328: Notes:
329: By default sets coloring type to IS_COLORING_GLOBAL
331: .seealso: MatColoringCreate(), ISColoringView(), ISColoringDestroy(), ISColoringSetType()
333: @*/
334: PetscErrorCode ISColoringCreate(MPI_Comm comm,PetscInt ncolors,PetscInt n,const ISColoringValue colors[],PetscCopyMode mode,ISColoring *iscoloring)
335: {
337: PetscMPIInt size,rank,tag;
338: PetscInt base,top,i;
339: PetscInt nc,ncwork;
340: MPI_Status status;
343: if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
344: 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);
345: 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);
346: }
347: PetscNew(iscoloring);
348: PetscCommDuplicate(comm,&(*iscoloring)->comm,&tag);
349: comm = (*iscoloring)->comm;
351: /* compute the number of the first node on my processor */
352: MPI_Comm_size(comm,&size);
354: /* should use MPI_Scan() */
355: MPI_Comm_rank(comm,&rank);
356: if (!rank) {
357: base = 0;
358: top = n;
359: } else {
360: MPI_Recv(&base,1,MPIU_INT,rank-1,tag,comm,&status);
361: top = base+n;
362: }
363: if (rank < size-1) {
364: MPI_Send(&top,1,MPIU_INT,rank+1,tag,comm);
365: }
367: /* compute the total number of colors */
368: ncwork = 0;
369: for (i=0; i<n; i++) {
370: if (ncwork < colors[i]) ncwork = colors[i];
371: }
372: ncwork++;
373: MPIU_Allreduce(&ncwork,&nc,1,MPIU_INT,MPI_MAX,comm);
374: 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);
375: (*iscoloring)->n = nc;
376: (*iscoloring)->is = 0;
377: (*iscoloring)->N = n;
378: (*iscoloring)->refct = 1;
379: (*iscoloring)->ctype = IS_COLORING_GLOBAL;
380: if (mode == PETSC_COPY_VALUES) {
381: PetscMalloc1(n,&(*iscoloring)->colors);
382: PetscLogObjectMemory((PetscObject)(*iscoloring),n*sizeof(ISColoringValue));
383: PetscArraycpy((*iscoloring)->colors,colors,n);
384: (*iscoloring)->allocated = PETSC_TRUE;
385: } else if (mode == PETSC_OWN_POINTER) {
386: (*iscoloring)->colors = (ISColoringValue*)colors;
387: (*iscoloring)->allocated = PETSC_TRUE;
388: } else {
389: (*iscoloring)->colors = (ISColoringValue*)colors;
390: (*iscoloring)->allocated = PETSC_FALSE;
391: }
392: ISColoringViewFromOptions(*iscoloring,NULL,"-is_coloring_view");
393: PetscInfo1(0,"Number of colors %D\n",nc);
394: return(0);
395: }
397: /*@
398: ISBuildTwoSided - Takes an IS that describes where we will go. Generates an IS that contains new numbers from remote or local
399: on the IS.
401: Collective on IS
403: Input Parameters
404: + ito - an IS describes where we will go. Negative target rank will be ignored
405: - toindx - an IS describes what indices should send. NULL means sending natural numbering
407: Output Parameter:
408: . rows - contains new numbers from remote or local
410: Level: advanced
412: .seealso: MatPartitioningCreate(), ISPartitioningToNumbering(), ISPartitioningCount()
414: @*/
415: PetscErrorCode ISBuildTwoSided(IS ito,IS toindx, IS *rows)
416: {
417: const PetscInt *ito_indices,*toindx_indices;
418: PetscInt *send_indices,rstart,*recv_indices,nrecvs,nsends;
419: PetscInt *tosizes,*fromsizes,i,j,*tosizes_tmp,*tooffsets_tmp,ito_ln;
420: PetscMPIInt *toranks,*fromranks,size,target_rank,*fromperm_newtoold,nto,nfrom;
421: PetscLayout isrmap;
422: MPI_Comm comm;
423: PetscSF sf;
424: PetscSFNode *iremote;
425: PetscErrorCode ierr;
428: PetscObjectGetComm((PetscObject)ito,&comm);
429: MPI_Comm_size(comm,&size);
430: ISGetLocalSize(ito,&ito_ln);
431: /* why we do not have ISGetLayout? */
432: isrmap = ito->map;
433: PetscLayoutGetRange(isrmap,&rstart,NULL);
434: ISGetIndices(ito,&ito_indices);
435: PetscCalloc2(size,&tosizes_tmp,size+1,&tooffsets_tmp);
436: for (i=0; i<ito_ln; i++) {
437: if (ito_indices[i]<0) continue;
438: #if defined(PETSC_USE_DEBUG)
439: if (ito_indices[i]>=size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"target rank %d is larger than communicator size %d ",ito_indices[i],size);
440: #endif
441: tosizes_tmp[ito_indices[i]]++;
442: }
443: nto = 0;
444: for (i=0; i<size; i++) {
445: tooffsets_tmp[i+1] = tooffsets_tmp[i]+tosizes_tmp[i];
446: if (tosizes_tmp[i]>0) nto++;
447: }
448: PetscCalloc2(nto,&toranks,2*nto,&tosizes);
449: nto = 0;
450: for (i=0; i<size; i++) {
451: if (tosizes_tmp[i]>0) {
452: toranks[nto] = i;
453: tosizes[2*nto] = tosizes_tmp[i];/* size */
454: tosizes[2*nto+1] = tooffsets_tmp[i];/* offset */
455: nto++;
456: }
457: }
458: nsends = tooffsets_tmp[size];
459: PetscCalloc1(nsends,&send_indices);
460: if (toindx) {
461: ISGetIndices(toindx,&toindx_indices);
462: }
463: for (i=0; i<ito_ln; i++) {
464: if (ito_indices[i]<0) continue;
465: target_rank = ito_indices[i];
466: send_indices[tooffsets_tmp[target_rank]] = toindx? toindx_indices[i]:(i+rstart);
467: tooffsets_tmp[target_rank]++;
468: }
469: if (toindx) {
470: ISRestoreIndices(toindx,&toindx_indices);
471: }
472: ISRestoreIndices(ito,&ito_indices);
473: PetscFree2(tosizes_tmp,tooffsets_tmp);
474: PetscCommBuildTwoSided(comm,2,MPIU_INT,nto,toranks,tosizes,&nfrom,&fromranks,&fromsizes);
475: PetscFree2(toranks,tosizes);
476: PetscMalloc1(nfrom,&fromperm_newtoold);
477: for (i=0; i<nfrom; i++) fromperm_newtoold[i] = i;
478: PetscSortMPIIntWithArray(nfrom,fromranks,fromperm_newtoold);
479: nrecvs = 0;
480: for (i=0; i<nfrom; i++) nrecvs += fromsizes[i*2];
481: PetscCalloc1(nrecvs,&recv_indices);
482: PetscMalloc1(nrecvs,&iremote);
483: nrecvs = 0;
484: for (i=0; i<nfrom; i++) {
485: for (j=0; j<fromsizes[2*fromperm_newtoold[i]]; j++) {
486: iremote[nrecvs].rank = fromranks[i];
487: iremote[nrecvs++].index = fromsizes[2*fromperm_newtoold[i]+1]+j;
488: }
489: }
490: PetscSFCreate(comm,&sf);
491: PetscSFSetGraph(sf,nsends,nrecvs,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
492: PetscSFSetType(sf,PETSCSFBASIC);
493: /* how to put a prefix ? */
494: PetscSFSetFromOptions(sf);
495: PetscSFBcastBegin(sf,MPIU_INT,send_indices,recv_indices);
496: PetscSFBcastEnd(sf,MPIU_INT,send_indices,recv_indices);
497: PetscSFDestroy(&sf);
498: PetscFree(fromranks);
499: PetscFree(fromsizes);
500: PetscFree(fromperm_newtoold);
501: PetscFree(send_indices);
502: if (rows) {
503: PetscSortInt(nrecvs,recv_indices);
504: ISCreateGeneral(comm,nrecvs,recv_indices,PETSC_OWN_POINTER,rows);
505: } else {
506: PetscFree(recv_indices);
507: }
508: return(0);
509: }
512: /*@
513: ISPartitioningToNumbering - Takes an ISPartitioning and on each processor
514: generates an IS that contains a new global node number for each index based
515: on the partitioing.
517: Collective on IS
519: Input Parameters
520: . partitioning - a partitioning as generated by MatPartitioningApply()
521: or MatPartitioningApplyND()
523: Output Parameter:
524: . is - on each processor the index set that defines the global numbers
525: (in the new numbering) for all the nodes currently (before the partitioning)
526: on that processor
528: Level: advanced
530: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningCount()
532: @*/
533: PetscErrorCode ISPartitioningToNumbering(IS part,IS *is)
534: {
535: MPI_Comm comm;
536: IS ndorder;
537: PetscInt i,np,npt,n,*starts = NULL,*sums = NULL,*lsizes = NULL,*newi = NULL;
538: const PetscInt *indices = NULL;
544: /* see if the partitioning comes from nested dissection */
545: PetscObjectQuery((PetscObject)part,"_petsc_matpartitioning_ndorder",(PetscObject*)&ndorder);
546: if (ndorder) {
547: PetscObjectReference((PetscObject)ndorder);
548: *is = ndorder;
549: return(0);
550: }
552: PetscObjectGetComm((PetscObject)part,&comm);
553: /* count the number of partitions, i.e., virtual processors */
554: ISGetLocalSize(part,&n);
555: ISGetIndices(part,&indices);
556: np = 0;
557: for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
558: MPIU_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
559: np = npt+1; /* so that it looks like a MPI_Comm_size output */
561: /*
562: lsizes - number of elements of each partition on this particular processor
563: sums - total number of "previous" nodes for any particular partition
564: starts - global number of first element in each partition on this processor
565: */
566: PetscMalloc3(np,&lsizes,np,&starts,np,&sums);
567: PetscArrayzero(lsizes,np);
568: for (i=0; i<n; i++) lsizes[indices[i]]++;
569: MPIU_Allreduce(lsizes,sums,np,MPIU_INT,MPI_SUM,comm);
570: MPI_Scan(lsizes,starts,np,MPIU_INT,MPI_SUM,comm);
571: for (i=0; i<np; i++) starts[i] -= lsizes[i];
572: for (i=1; i<np; i++) {
573: sums[i] += sums[i-1];
574: starts[i] += sums[i-1];
575: }
577: /*
578: For each local index give it the new global number
579: */
580: PetscMalloc1(n,&newi);
581: for (i=0; i<n; i++) newi[i] = starts[indices[i]]++;
582: PetscFree3(lsizes,starts,sums);
584: ISRestoreIndices(part,&indices);
585: ISCreateGeneral(comm,n,newi,PETSC_OWN_POINTER,is);
586: ISSetPermutation(*is);
587: return(0);
588: }
590: /*@
591: ISPartitioningCount - Takes a ISPartitioning and determines the number of
592: resulting elements on each (partition) process
594: Collective on IS
596: Input Parameters:
597: + partitioning - a partitioning as generated by MatPartitioningApply() or
598: MatPartitioningApplyND()
599: - len - length of the array count, this is the total number of partitions
601: Output Parameter:
602: . count - array of length size, to contain the number of elements assigned
603: to each partition, where size is the number of partitions generated
604: (see notes below).
606: Level: advanced
608: Notes:
609: By default the number of partitions generated (and thus the length
610: of count) is the size of the communicator associated with IS,
611: but it can be set by MatPartitioningSetNParts. The resulting array
612: of lengths can for instance serve as input of PCBJacobiSetTotalBlocks.
613: If the partitioning has been obtained by MatPartitioningApplyND(),
614: the returned count does not include the separators.
616: .seealso: MatPartitioningCreate(), AOCreateBasic(), ISPartitioningToNumbering(),
617: MatPartitioningSetNParts(), MatPartitioningApply(), MatPartitioningApplyND()
619: @*/
620: PetscErrorCode ISPartitioningCount(IS part,PetscInt len,PetscInt count[])
621: {
622: MPI_Comm comm;
623: PetscInt i,n,*lsizes;
624: const PetscInt *indices;
626: PetscMPIInt npp;
629: PetscObjectGetComm((PetscObject)part,&comm);
630: if (len == PETSC_DEFAULT) {
631: PetscMPIInt size;
632: MPI_Comm_size(comm,&size);
633: len = (PetscInt) size;
634: }
636: /* count the number of partitions */
637: ISGetLocalSize(part,&n);
638: ISGetIndices(part,&indices);
639: #if defined(PETSC_USE_DEBUG)
640: {
641: PetscInt np = 0,npt;
642: for (i=0; i<n; i++) np = PetscMax(np,indices[i]);
643: MPIU_Allreduce(&np,&npt,1,MPIU_INT,MPI_MAX,comm);
644: np = npt+1; /* so that it looks like a MPI_Comm_size output */
645: 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);
646: }
647: #endif
649: /*
650: lsizes - number of elements of each partition on this particular processor
651: sums - total number of "previous" nodes for any particular partition
652: starts - global number of first element in each partition on this processor
653: */
654: PetscCalloc1(len,&lsizes);
655: for (i=0; i<n; i++) {
656: if (indices[i] > -1) lsizes[indices[i]]++;
657: }
658: ISRestoreIndices(part,&indices);
659: PetscMPIIntCast(len,&npp);
660: MPIU_Allreduce(lsizes,count,npp,MPIU_INT,MPI_SUM,comm);
661: PetscFree(lsizes);
662: return(0);
663: }
665: /*@
666: ISAllGather - Given an index set (IS) on each processor, generates a large
667: index set (same on each processor) by concatenating together each
668: processors index set.
670: Collective on IS
672: Input Parameter:
673: . is - the distributed index set
675: Output Parameter:
676: . isout - the concatenated index set (same on all processors)
678: Notes:
679: ISAllGather() is clearly not scalable for large index sets.
681: The IS created on each processor must be created with a common
682: communicator (e.g., PETSC_COMM_WORLD). If the index sets were created
683: with PETSC_COMM_SELF, this routine will not work as expected, since
684: each process will generate its own new IS that consists only of
685: itself.
687: The communicator for this new IS is PETSC_COMM_SELF
689: Level: intermediate
691: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock()
692: @*/
693: PetscErrorCode ISAllGather(IS is,IS *isout)
694: {
696: PetscInt *indices,n,i,N,step,first;
697: const PetscInt *lindices;
698: MPI_Comm comm;
699: PetscMPIInt size,*sizes = NULL,*offsets = NULL,nn;
700: PetscBool stride;
706: PetscObjectGetComm((PetscObject)is,&comm);
707: MPI_Comm_size(comm,&size);
708: ISGetLocalSize(is,&n);
709: PetscObjectTypeCompare((PetscObject)is,ISSTRIDE,&stride);
710: if (size == 1 && stride) { /* should handle parallel ISStride also */
711: ISStrideGetInfo(is,&first,&step);
712: ISCreateStride(PETSC_COMM_SELF,n,first,step,isout);
713: } else {
714: PetscMalloc2(size,&sizes,size,&offsets);
716: PetscMPIIntCast(n,&nn);
717: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
718: offsets[0] = 0;
719: for (i=1; i<size; i++) {
720: PetscInt s = offsets[i-1] + sizes[i-1];
721: PetscMPIIntCast(s,&offsets[i]);
722: }
723: N = offsets[size-1] + sizes[size-1];
725: PetscMalloc1(N,&indices);
726: ISGetIndices(is,&lindices);
727: MPI_Allgatherv((void*)lindices,nn,MPIU_INT,indices,sizes,offsets,MPIU_INT,comm);
728: ISRestoreIndices(is,&lindices);
729: PetscFree2(sizes,offsets);
731: ISCreateGeneral(PETSC_COMM_SELF,N,indices,PETSC_OWN_POINTER,isout);
732: }
733: return(0);
734: }
736: /*@C
737: ISAllGatherColors - Given a a set of colors on each processor, generates a large
738: set (same on each processor) by concatenating together each processors colors
740: Collective
742: Input Parameter:
743: + comm - communicator to share the indices
744: . n - local size of set
745: - lindices - local colors
747: Output Parameter:
748: + outN - total number of indices
749: - outindices - all of the colors
751: Notes:
752: ISAllGatherColors() is clearly not scalable for large index sets.
755: Level: intermediate
757: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
758: @*/
759: PetscErrorCode ISAllGatherColors(MPI_Comm comm,PetscInt n,ISColoringValue *lindices,PetscInt *outN,ISColoringValue *outindices[])
760: {
761: ISColoringValue *indices;
762: PetscErrorCode ierr;
763: PetscInt i,N;
764: PetscMPIInt size,*offsets = NULL,*sizes = NULL, nn = n;
767: MPI_Comm_size(comm,&size);
768: PetscMalloc2(size,&sizes,size,&offsets);
770: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
771: offsets[0] = 0;
772: for (i=1; i<size; i++) offsets[i] = offsets[i-1] + sizes[i-1];
773: N = offsets[size-1] + sizes[size-1];
774: PetscFree2(sizes,offsets);
776: PetscMalloc1(N+1,&indices);
777: MPI_Allgatherv(lindices,(PetscMPIInt)n,MPIU_COLORING_VALUE,indices,sizes,offsets,MPIU_COLORING_VALUE,comm);
779: *outindices = indices;
780: if (outN) *outN = N;
781: return(0);
782: }
784: /*@
785: ISComplement - Given an index set (IS) generates the complement index set. That is all
786: all indices that are NOT in the given set.
788: Collective on IS
790: Input Parameter:
791: + is - the index set
792: . nmin - the first index desired in the local part of the complement
793: - 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)
795: Output Parameter:
796: . isout - the complement
798: Notes:
799: The communicator for this new IS is the same as for the input IS
801: For a parallel IS, this will generate the local part of the complement on each process
803: To generate the entire complement (on each process) of a parallel IS, first call ISAllGather() and then
804: call this routine.
806: Level: intermediate
808: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
809: @*/
810: PetscErrorCode ISComplement(IS is,PetscInt nmin,PetscInt nmax,IS *isout)
811: {
813: const PetscInt *indices;
814: PetscInt n,i,j,unique,cnt,*nindices;
815: PetscBool sorted;
820: if (nmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nmin %D cannot be negative",nmin);
821: if (nmin > nmax) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nmin %D cannot be greater than nmax %D",nmin,nmax);
822: ISSorted(is,&sorted);
823: if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Index set must be sorted");
825: ISGetLocalSize(is,&n);
826: ISGetIndices(is,&indices);
827: #if defined(PETSC_USE_DEBUG)
828: for (i=0; i<n; i++) {
829: 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);
830: 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);
831: }
832: #endif
833: /* Count number of unique entries */
834: unique = (n>0);
835: for (i=0; i<n-1; i++) {
836: if (indices[i+1] != indices[i]) unique++;
837: }
838: PetscMalloc1(nmax-nmin-unique,&nindices);
839: cnt = 0;
840: for (i=nmin,j=0; i<nmax; i++) {
841: if (j<n && i==indices[j]) do { j++; } while (j<n && i==indices[j]);
842: else nindices[cnt++] = i;
843: }
844: 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);
845: ISCreateGeneral(PetscObjectComm((PetscObject)is),cnt,nindices,PETSC_OWN_POINTER,isout);
846: ISRestoreIndices(is,&indices);
847: return(0);
848: }