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: }