1 
2 /*
3      Provides utility routines for split MPI communicator.
4 */
5 #include <petscsys.h>    /*I   "petscsys.h"    I*/
6 #include <petscviewer.h>
7 
8 const char *const PetscSubcommTypes[] = {"GENERAL","CONTIGUOUS","INTERLACED","PetscSubcommType","PETSC_SUBCOMM_",NULL};
9 
10 static PetscErrorCode PetscSubcommCreate_contiguous(PetscSubcomm);
11 static PetscErrorCode PetscSubcommCreate_interlaced(PetscSubcomm);
12 
13 /*@
14    PetscSubcommSetFromOptions - Allows setting options from a PetscSubcomm
15 
16    Collective on PetscSubcomm
17 
18    Input Parameter:
19 .  psubcomm - PetscSubcomm context
20 
21    Level: beginner
22 
23 @*/
PetscSubcommSetFromOptions(PetscSubcomm psubcomm)24 PetscErrorCode PetscSubcommSetFromOptions(PetscSubcomm psubcomm)
25 {
26   PetscErrorCode   ierr;
27   PetscSubcommType type;
28   PetscBool        flg;
29 
30   PetscFunctionBegin;
31   if (!psubcomm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Must call PetscSubcommCreate firt");
32 
33   ierr = PetscOptionsBegin(psubcomm->parent,psubcomm->subcommprefix,"Options for PetscSubcomm",NULL);CHKERRQ(ierr);
34   ierr = PetscOptionsEnum("-psubcomm_type",NULL,NULL,PetscSubcommTypes,(PetscEnum)psubcomm->type,(PetscEnum*)&type,&flg);CHKERRQ(ierr);
35   if (flg && psubcomm->type != type) {
36     /* free old structures */
37     ierr = PetscCommDestroy(&(psubcomm)->dupparent);CHKERRQ(ierr);
38     ierr = PetscCommDestroy(&(psubcomm)->child);CHKERRQ(ierr);
39     ierr = PetscFree((psubcomm)->subsize);CHKERRQ(ierr);
40     switch (type) {
41     case PETSC_SUBCOMM_GENERAL:
42       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Runtime option PETSC_SUBCOMM_GENERAL is not supported, use PetscSubcommSetTypeGeneral()");
43     case PETSC_SUBCOMM_CONTIGUOUS:
44       ierr = PetscSubcommCreate_contiguous(psubcomm);CHKERRQ(ierr);
45       break;
46     case PETSC_SUBCOMM_INTERLACED:
47       ierr = PetscSubcommCreate_interlaced(psubcomm);CHKERRQ(ierr);
48       break;
49     default:
50       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"PetscSubcommType %s is not supported yet",PetscSubcommTypes[type]);
51     }
52   }
53 
54   ierr = PetscOptionsName("-psubcomm_view","Triggers display of PetscSubcomm context","PetscSubcommView",&flg);CHKERRQ(ierr);
55   if (flg) {
56     ierr = PetscSubcommView(psubcomm,PETSC_VIEWER_STDOUT_(psubcomm->parent));CHKERRQ(ierr);
57   }
58   ierr = PetscOptionsEnd();CHKERRQ(ierr);
59   PetscFunctionReturn(0);
60 }
61 
62 /*@C
63   PetscSubcommSetOptionsPrefix - Sets the prefix used for searching for all
64   PetscSubcomm items in the options database.
65 
66   Logically collective on PetscSubcomm.
67 
68   Level: Intermediate
69 
70   Input Parameters:
71 +   psubcomm - PetscSubcomm context
72 -   prefix - the prefix to prepend all PetscSubcomm item names with.
73 
74 @*/
PetscSubcommSetOptionsPrefix(PetscSubcomm psubcomm,const char pre[])75 PetscErrorCode PetscSubcommSetOptionsPrefix(PetscSubcomm psubcomm,const char pre[])
76 {
77   PetscErrorCode   ierr;
78 
79   PetscFunctionBegin;
80    if (!pre) {
81     ierr = PetscFree(psubcomm->subcommprefix);CHKERRQ(ierr);
82   } else {
83     if (pre[0] == '-') SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Options prefix should not begin with a hypen");
84     ierr = PetscFree(psubcomm->subcommprefix);CHKERRQ(ierr);
85     ierr = PetscStrallocpy(pre,&(psubcomm->subcommprefix));CHKERRQ(ierr);
86   }
87   PetscFunctionReturn(0);
88 }
89 
90 /*@C
91    PetscSubcommView - Views a PetscSubcomm of values as either ASCII text or a binary file
92 
93    Collective on PetscSubcomm
94 
95    Input Parameter:
96 +  psubcomm - PetscSubcomm context
97 -  viewer - location to view the values
98 
99    Level: beginner
100 @*/
PetscSubcommView(PetscSubcomm psubcomm,PetscViewer viewer)101 PetscErrorCode PetscSubcommView(PetscSubcomm psubcomm,PetscViewer viewer)
102 {
103   PetscErrorCode    ierr;
104   PetscBool         iascii;
105   PetscViewerFormat format;
106 
107   PetscFunctionBegin;
108   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
109   if (iascii) {
110     ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
111     if (format == PETSC_VIEWER_DEFAULT) {
112       MPI_Comm    comm=psubcomm->parent;
113       PetscMPIInt rank,size,subsize,subrank,duprank;
114 
115       ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
116       ierr = PetscViewerASCIIPrintf(viewer,"PetscSubcomm type %s with total %d MPI processes:\n",PetscSubcommTypes[psubcomm->type],size);CHKERRQ(ierr);
117       ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
118       ierr = MPI_Comm_size(psubcomm->child,&subsize);CHKERRQ(ierr);
119       ierr = MPI_Comm_rank(psubcomm->child,&subrank);CHKERRQ(ierr);
120       ierr = MPI_Comm_rank(psubcomm->dupparent,&duprank);CHKERRQ(ierr);
121       ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr);
122       ierr = PetscViewerASCIISynchronizedPrintf(viewer,"  [%d], color %d, sub-size %d, sub-rank %d, duprank %d\n",rank,psubcomm->color,subsize,subrank,duprank);CHKERRQ(ierr);
123       ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
124       ierr = PetscViewerASCIIPopSynchronized(viewer);CHKERRQ(ierr);
125     }
126   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not supported yet");
127   PetscFunctionReturn(0);
128 }
129 
130 /*@
131   PetscSubcommSetNumber - Set total number of subcommunicators.
132 
133    Collective
134 
135    Input Parameter:
136 +  psubcomm - PetscSubcomm context
137 -  nsubcomm - the total number of subcommunicators in psubcomm
138 
139    Level: advanced
140 
141 .seealso: PetscSubcommCreate(),PetscSubcommDestroy(),PetscSubcommSetType(),PetscSubcommSetTypeGeneral()
142 @*/
PetscSubcommSetNumber(PetscSubcomm psubcomm,PetscInt nsubcomm)143 PetscErrorCode  PetscSubcommSetNumber(PetscSubcomm psubcomm,PetscInt nsubcomm)
144 {
145   PetscErrorCode ierr;
146   MPI_Comm       comm=psubcomm->parent;
147   PetscMPIInt    msub,size;
148 
149   PetscFunctionBegin;
150   if (!psubcomm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"PetscSubcomm is not created. Call PetscSubcommCreate() first");
151   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
152   ierr = PetscMPIIntCast(nsubcomm,&msub);CHKERRQ(ierr);
153   if (msub < 1 || msub > size) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Num of subcommunicators %d cannot be < 1 or > input comm size %d",msub,size);
154 
155   psubcomm->n = msub;
156   PetscFunctionReturn(0);
157 }
158 
159 /*@
160   PetscSubcommSetType - Set type of subcommunicators.
161 
162    Collective
163 
164    Input Parameter:
165 +  psubcomm - PetscSubcomm context
166 -  subcommtype - subcommunicator type, PETSC_SUBCOMM_CONTIGUOUS,PETSC_SUBCOMM_INTERLACED
167 
168    Level: advanced
169 
170 .seealso: PetscSubcommCreate(),PetscSubcommDestroy(),PetscSubcommSetNumber(),PetscSubcommSetTypeGeneral(), PetscSubcommType
171 @*/
PetscSubcommSetType(PetscSubcomm psubcomm,PetscSubcommType subcommtype)172 PetscErrorCode  PetscSubcommSetType(PetscSubcomm psubcomm,PetscSubcommType subcommtype)
173 {
174   PetscErrorCode ierr;
175 
176   PetscFunctionBegin;
177   if (!psubcomm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"PetscSubcomm is not created. Call PetscSubcommCreate()");
178   if (psubcomm->n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of subcommunicators %d is incorrect. Call PetscSubcommSetNumber()",psubcomm->n);
179 
180   if (subcommtype == PETSC_SUBCOMM_CONTIGUOUS) {
181     ierr = PetscSubcommCreate_contiguous(psubcomm);CHKERRQ(ierr);
182   } else if (subcommtype == PETSC_SUBCOMM_INTERLACED) {
183     ierr = PetscSubcommCreate_interlaced(psubcomm);CHKERRQ(ierr);
184   } else SETERRQ1(psubcomm->parent,PETSC_ERR_SUP,"PetscSubcommType %s is not supported yet",PetscSubcommTypes[subcommtype]);
185   PetscFunctionReturn(0);
186 }
187 
188 /*@
189   PetscSubcommSetTypeGeneral - Set a PetscSubcomm from user's specifications
190 
191    Collective
192 
193    Input Parameter:
194 +  psubcomm - PetscSubcomm context
195 .  color   - control of subset assignment (nonnegative integer). Processes with the same color are in the same subcommunicator.
196 -  subrank - rank in the subcommunicator
197 
198    Level: advanced
199 
200 .seealso: PetscSubcommCreate(),PetscSubcommDestroy(),PetscSubcommSetNumber(),PetscSubcommSetType()
201 @*/
PetscSubcommSetTypeGeneral(PetscSubcomm psubcomm,PetscMPIInt color,PetscMPIInt subrank)202 PetscErrorCode PetscSubcommSetTypeGeneral(PetscSubcomm psubcomm,PetscMPIInt color,PetscMPIInt subrank)
203 {
204   PetscErrorCode ierr;
205   MPI_Comm       subcomm=0,dupcomm=0,comm=psubcomm->parent;
206   PetscMPIInt    size,icolor,duprank,*recvbuf,sendbuf[3],mysubsize,rank,*subsize;
207   PetscMPIInt    i,nsubcomm=psubcomm->n;
208 
209   PetscFunctionBegin;
210   if (!psubcomm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"PetscSubcomm is not created. Call PetscSubcommCreate()");
211   if (nsubcomm < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of subcommunicators %d is incorrect. Call PetscSubcommSetNumber()",nsubcomm);
212 
213   ierr = MPI_Comm_split(comm,color,subrank,&subcomm);CHKERRQ(ierr);
214 
215   /* create dupcomm with same size as comm, but its rank, duprank, maps subcomm's contiguously into dupcomm */
216   /* TODO: this can be done in an ostensibly scalale way (i.e., without allocating an array of size 'size') as is done in PetscObjectsCreateGlobalOrdering(). */
217   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
218   ierr = PetscMalloc1(2*size,&recvbuf);CHKERRQ(ierr);
219 
220   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
221   ierr = MPI_Comm_size(subcomm,&mysubsize);CHKERRQ(ierr);
222 
223   sendbuf[0] = color;
224   sendbuf[1] = mysubsize;
225   ierr = MPI_Allgather(sendbuf,2,MPI_INT,recvbuf,2,MPI_INT,comm);CHKERRQ(ierr);
226 
227   ierr = PetscCalloc1(nsubcomm,&subsize);CHKERRQ(ierr);
228   for (i=0; i<2*size; i+=2) {
229     subsize[recvbuf[i]] = recvbuf[i+1];
230   }
231   ierr = PetscFree(recvbuf);CHKERRQ(ierr);
232 
233   duprank = 0;
234   for (icolor=0; icolor<nsubcomm; icolor++) {
235     if (icolor != color) { /* not color of this process */
236       duprank += subsize[icolor];
237     } else {
238       duprank += subrank;
239       break;
240     }
241   }
242   ierr = MPI_Comm_split(comm,0,duprank,&dupcomm);CHKERRQ(ierr);
243 
244   ierr = PetscCommDuplicate(dupcomm,&psubcomm->dupparent,NULL);CHKERRQ(ierr);
245   ierr = PetscCommDuplicate(subcomm,&psubcomm->child,NULL);CHKERRQ(ierr);
246   ierr = MPI_Comm_free(&dupcomm);CHKERRQ(ierr);
247   ierr = MPI_Comm_free(&subcomm);CHKERRQ(ierr);
248 
249   psubcomm->color   = color;
250   psubcomm->subsize = subsize;
251   psubcomm->type    = PETSC_SUBCOMM_GENERAL;
252   PetscFunctionReturn(0);
253 }
254 
255 /*@
256   PetscSubcommDestroy - Destroys a PetscSubcomm object
257 
258    Collective on PetscSubcomm
259 
260    Input Parameter:
261    .  psubcomm - the PetscSubcomm context
262 
263    Level: advanced
264 
265 .seealso: PetscSubcommCreate(),PetscSubcommSetType()
266 @*/
PetscSubcommDestroy(PetscSubcomm * psubcomm)267 PetscErrorCode  PetscSubcommDestroy(PetscSubcomm *psubcomm)
268 {
269   PetscErrorCode ierr;
270 
271   PetscFunctionBegin;
272   if (!*psubcomm) PetscFunctionReturn(0);
273   ierr = PetscCommDestroy(&(*psubcomm)->dupparent);CHKERRQ(ierr);
274   ierr = PetscCommDestroy(&(*psubcomm)->child);CHKERRQ(ierr);
275   ierr = PetscFree((*psubcomm)->subsize);CHKERRQ(ierr);
276   if ((*psubcomm)->subcommprefix) { ierr = PetscFree((*psubcomm)->subcommprefix);CHKERRQ(ierr); }
277   ierr = PetscFree((*psubcomm));CHKERRQ(ierr);
278   PetscFunctionReturn(0);
279 }
280 
281 /*@
282   PetscSubcommCreate - Create a PetscSubcomm context.
283 
284    Collective
285 
286    Input Parameter:
287 .  comm - MPI communicator
288 
289    Output Parameter:
290 .  psubcomm - location to store the PetscSubcomm context
291 
292    Level: advanced
293 
294 .seealso: PetscSubcommDestroy(), PetscSubcommSetTypeGeneral(), PetscSubcommSetFromOptions(), PetscSubcommSetType(),
295           PetscSubcommSetNumber()
296 @*/
PetscSubcommCreate(MPI_Comm comm,PetscSubcomm * psubcomm)297 PetscErrorCode  PetscSubcommCreate(MPI_Comm comm,PetscSubcomm *psubcomm)
298 {
299   PetscErrorCode ierr;
300   PetscMPIInt    rank,size;
301 
302   PetscFunctionBegin;
303   ierr = PetscNew(psubcomm);CHKERRQ(ierr);
304 
305   /* set defaults */
306   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
307   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
308 
309   (*psubcomm)->parent    = comm;
310   (*psubcomm)->dupparent = comm;
311   (*psubcomm)->child     = PETSC_COMM_SELF;
312   (*psubcomm)->n         = size;
313   (*psubcomm)->color     = rank;
314   (*psubcomm)->subsize   = NULL;
315   (*psubcomm)->type      = PETSC_SUBCOMM_INTERLACED;
316   PetscFunctionReturn(0);
317 }
318 
319 /*@C
320   PetscSubcommGetParent - Gets the communicator that was used to create the PetscSubcomm
321 
322    Collective
323 
324    Input Parameter:
325 .  scomm - the PetscSubcomm
326 
327    Output Parameter:
328 .  pcomm - location to store the parent communicator
329 
330    Level: intermediate
331 
332 .seealso: PetscSubcommDestroy(), PetscSubcommSetTypeGeneral(), PetscSubcommSetFromOptions(), PetscSubcommSetType(),
333           PetscSubcommSetNumber(), PetscSubcommGetChild(), PetscSubcommContiguousParent()
334 @*/
PetscSubcommGetParent(PetscSubcomm scomm,MPI_Comm * pcomm)335 PetscErrorCode  PetscSubcommGetParent(PetscSubcomm scomm,MPI_Comm *pcomm)
336 {
337   *pcomm = PetscSubcommParent(scomm);
338   return 0;
339 }
340 
341 /*@C
342   PetscSubcommGetContiguousParent - Gets a communicator that that is a duplicate of the parent but has the ranks
343                                     reordered by the order they are in the children
344 
345    Collective
346 
347    Input Parameter:
348 .  scomm - the PetscSubcomm
349 
350    Output Parameter:
351 .  pcomm - location to store the parent communicator
352 
353    Level: intermediate
354 
355 .seealso: PetscSubcommDestroy(), PetscSubcommSetTypeGeneral(), PetscSubcommSetFromOptions(), PetscSubcommSetType(),
356           PetscSubcommSetNumber(), PetscSubcommGetChild(), PetscSubcommContiguousParent()
357 @*/
PetscSubcommGetContiguousParent(PetscSubcomm scomm,MPI_Comm * pcomm)358 PetscErrorCode  PetscSubcommGetContiguousParent(PetscSubcomm scomm,MPI_Comm *pcomm)
359 {
360   *pcomm = PetscSubcommContiguousParent(scomm);
361   return 0;
362 }
363 
364 /*@C
365   PetscSubcommGetChild - Gets the communicator created by the PetscSubcomm
366 
367    Collective
368 
369    Input Parameter:
370 .  scomm - the PetscSubcomm
371 
372    Output Parameter:
373 .  ccomm - location to store the child communicator
374 
375    Level: intermediate
376 
377 .seealso: PetscSubcommDestroy(), PetscSubcommSetTypeGeneral(), PetscSubcommSetFromOptions(), PetscSubcommSetType(),
378           PetscSubcommSetNumber(), PetscSubcommGetParent(), PetscSubcommContiguousParent()
379 @*/
PetscSubcommGetChild(PetscSubcomm scomm,MPI_Comm * ccomm)380 PetscErrorCode  PetscSubcommGetChild(PetscSubcomm scomm,MPI_Comm *ccomm)
381 {
382   *ccomm = PetscSubcommChild(scomm);
383   return 0;
384 }
385 
PetscSubcommCreate_contiguous(PetscSubcomm psubcomm)386 static PetscErrorCode PetscSubcommCreate_contiguous(PetscSubcomm psubcomm)
387 {
388   PetscErrorCode ierr;
389   PetscMPIInt    rank,size,*subsize,duprank=-1,subrank=-1;
390   PetscMPIInt    np_subcomm,nleftover,i,color=-1,rankstart,nsubcomm=psubcomm->n;
391   MPI_Comm       subcomm=0,dupcomm=0,comm=psubcomm->parent;
392 
393   PetscFunctionBegin;
394   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
395   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
396 
397   /* get size of each subcommunicator */
398   ierr = PetscMalloc1(1+nsubcomm,&subsize);CHKERRQ(ierr);
399 
400   np_subcomm = size/nsubcomm;
401   nleftover  = size - nsubcomm*np_subcomm;
402   for (i=0; i<nsubcomm; i++) {
403     subsize[i] = np_subcomm;
404     if (i<nleftover) subsize[i]++;
405   }
406 
407   /* get color and subrank of this proc */
408   rankstart = 0;
409   for (i=0; i<nsubcomm; i++) {
410     if (rank >= rankstart && rank < rankstart+subsize[i]) {
411       color   = i;
412       subrank = rank - rankstart;
413       duprank = rank;
414       break;
415     } else rankstart += subsize[i];
416   }
417 
418   ierr = MPI_Comm_split(comm,color,subrank,&subcomm);CHKERRQ(ierr);
419 
420   /* create dupcomm with same size as comm, but its rank, duprank, maps subcomm's contiguously into dupcomm */
421   ierr = MPI_Comm_split(comm,0,duprank,&dupcomm);CHKERRQ(ierr);
422   ierr = PetscCommDuplicate(dupcomm,&psubcomm->dupparent,NULL);CHKERRQ(ierr);
423   ierr = PetscCommDuplicate(subcomm,&psubcomm->child,NULL);CHKERRQ(ierr);
424   ierr = MPI_Comm_free(&dupcomm);CHKERRQ(ierr);
425   ierr = MPI_Comm_free(&subcomm);CHKERRQ(ierr);
426 
427   psubcomm->color   = color;
428   psubcomm->subsize = subsize;
429   psubcomm->type    = PETSC_SUBCOMM_CONTIGUOUS;
430   PetscFunctionReturn(0);
431 }
432 
433 /*
434    Note:
435    In PCREDUNDANT, to avoid data scattering from subcomm back to original comm, we create subcommunicators
436    by iteratively taking a process into a subcommunicator.
437    Example: size=4, nsubcomm=(*psubcomm)->n=3
438      comm=(*psubcomm)->parent:
439       rank:     [0]  [1]  [2]  [3]
440       color:     0    1    2    0
441 
442      subcomm=(*psubcomm)->comm:
443       subrank:  [0]  [0]  [0]  [1]
444 
445      dupcomm=(*psubcomm)->dupparent:
446       duprank:  [0]  [2]  [3]  [1]
447 
448      Here, subcomm[color = 0] has subsize=2, owns process [0] and [3]
449            subcomm[color = 1] has subsize=1, owns process [1]
450            subcomm[color = 2] has subsize=1, owns process [2]
451            dupcomm has same number of processes as comm, and its duprank maps
452            processes in subcomm contiguously into a 1d array:
453             duprank: [0] [1]      [2]         [3]
454             rank:    [0] [3]      [1]         [2]
455                     subcomm[0] subcomm[1]  subcomm[2]
456 */
457 
PetscSubcommCreate_interlaced(PetscSubcomm psubcomm)458 static PetscErrorCode PetscSubcommCreate_interlaced(PetscSubcomm psubcomm)
459 {
460   PetscErrorCode ierr;
461   PetscMPIInt    rank,size,*subsize,duprank,subrank;
462   PetscMPIInt    np_subcomm,nleftover,i,j,color,nsubcomm=psubcomm->n;
463   MPI_Comm       subcomm=0,dupcomm=0,comm=psubcomm->parent;
464 
465   PetscFunctionBegin;
466   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
467   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
468 
469   /* get size of each subcommunicator */
470   ierr = PetscMalloc1(1+nsubcomm,&subsize);CHKERRQ(ierr);
471 
472   np_subcomm = size/nsubcomm;
473   nleftover  = size - nsubcomm*np_subcomm;
474   for (i=0; i<nsubcomm; i++) {
475     subsize[i] = np_subcomm;
476     if (i<nleftover) subsize[i]++;
477   }
478 
479   /* find color for this proc */
480   color   = rank%nsubcomm;
481   subrank = rank/nsubcomm;
482 
483   ierr = MPI_Comm_split(comm,color,subrank,&subcomm);CHKERRQ(ierr);
484 
485   j = 0; duprank = 0;
486   for (i=0; i<nsubcomm; i++) {
487     if (j == color) {
488       duprank += subrank;
489       break;
490     }
491     duprank += subsize[i]; j++;
492   }
493 
494   /* create dupcomm with same size as comm, but its rank, duprank, maps subcomm's contiguously into dupcomm */
495   ierr = MPI_Comm_split(comm,0,duprank,&dupcomm);CHKERRQ(ierr);
496   ierr = PetscCommDuplicate(dupcomm,&psubcomm->dupparent,NULL);CHKERRQ(ierr);
497   ierr = PetscCommDuplicate(subcomm,&psubcomm->child,NULL);CHKERRQ(ierr);
498   ierr = MPI_Comm_free(&dupcomm);CHKERRQ(ierr);
499   ierr = MPI_Comm_free(&subcomm);CHKERRQ(ierr);
500 
501   psubcomm->color   = color;
502   psubcomm->subsize = subsize;
503   psubcomm->type    = PETSC_SUBCOMM_INTERLACED;
504   PetscFunctionReturn(0);
505 }
506 
507 
508