Actual source code: isltog.c

  1: #define PETSCVEC_DLL

 3:  #include petscvec.h
 4:  #include src/vec/is/isimpl.h

  6: PetscCookie  IS_LTOGM_COOKIE = -1;

 10: /*@C
 11:     ISLocalToGlobalMappingGetSize - Gets the local size of a local to global mapping.

 13:     Not Collective

 15:     Input Parameter:
 16: .   ltog - local to global mapping

 18:     Output Parameter:
 19: .   n - the number of entries in the local mapping

 21:     Level: advanced

 23:     Concepts: mapping^local to global

 25: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate()
 26: @*/
 27: PetscErrorCode  ISLocalToGlobalMappingGetSize(ISLocalToGlobalMapping mapping,PetscInt *n)
 28: {
 32:   *n = mapping->n;
 33:   return(0);
 34: }

 38: /*@C
 39:     ISLocalToGlobalMappingView - View a local to global mapping

 41:     Not Collective

 43:     Input Parameters:
 44: +   ltog - local to global mapping
 45: -   viewer - viewer

 47:     Level: advanced

 49:     Concepts: mapping^local to global

 51: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate()
 52: @*/
 53: PetscErrorCode  ISLocalToGlobalMappingView(ISLocalToGlobalMapping mapping,PetscViewer viewer)
 54: {
 55:   PetscInt        i;
 56:   PetscMPIInt     rank;
 57:   PetscTruth      iascii;
 58:   PetscErrorCode  ierr;

 62:   if (!viewer) {
 63:     PetscViewerASCIIGetStdout(mapping->comm,&viewer);
 64:   }

 67:   MPI_Comm_rank(mapping->comm,&rank);
 68:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
 69:   if (iascii) {
 70:     for (i=0; i<mapping->n; i++) {
 71:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %d %d\n",rank,i,mapping->indices[i]);
 72:     }
 73:     PetscViewerFlush(viewer);
 74:   } else {
 75:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for ISLocalToGlobalMapping",((PetscObject)viewer)->type_name);
 76:   }

 78:   return(0);
 79: }

 83: /*@
 84:     ISLocalToGlobalMappingCreateIS - Creates a mapping between a local (0 to n)
 85:     ordering and a global parallel ordering.

 87:     Not collective

 89:     Input Parameter:
 90: .   is - index set containing the global numbers for each local

 92:     Output Parameter:
 93: .   mapping - new mapping data structure

 95:     Level: advanced

 97:     Concepts: mapping^local to global

 99: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate()
100: @*/
101: PetscErrorCode  ISLocalToGlobalMappingCreateIS(IS is,ISLocalToGlobalMapping *mapping)
102: {
104:   PetscInt      n,*indices;
105:   MPI_Comm comm;


111:   PetscObjectGetComm((PetscObject)is,&comm);
112:   ISGetLocalSize(is,&n);
113:   ISGetIndices(is,&indices);
114:   ISLocalToGlobalMappingCreate(comm,n,indices,mapping);
115:   ISRestoreIndices(is,&indices);

117:   return(0);
118: }


123: /*@
124:     ISLocalToGlobalMappingCreate - Creates a mapping between a local (0 to n)
125:     ordering and a global parallel ordering.

127:     Not Collective, but communicator may have more than one process

129:     Input Parameters:
130: +   comm - MPI communicator
131: .   n - the number of local elements
132: -   indices - the global index for each local element

134:     Output Parameter:
135: .   mapping - new mapping data structure

137:     Level: advanced

139:     Concepts: mapping^local to global

141: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreateNC()
142: @*/
143: PetscErrorCode  ISLocalToGlobalMappingCreate(MPI_Comm cm,PetscInt n,const PetscInt indices[],ISLocalToGlobalMapping *mapping)
144: {
146:   PetscInt       *in;

151:   PetscMalloc(n*sizeof(PetscInt),&in);
152:   PetscMemcpy(in,indices,n*sizeof(PetscInt));
153:   ISLocalToGlobalMappingCreateNC(cm,n,in,mapping);
154:   return(0);
155: }

159: /*@C
160:     ISLocalToGlobalMappingCreateNC - Creates a mapping between a local (0 to n)
161:     ordering and a global parallel ordering.

163:     Not Collective, but communicator may have more than one process

165:     Input Parameters:
166: +   comm - MPI communicator
167: .   n - the number of local elements
168: -   indices - the global index for each local element

170:     Output Parameter:
171: .   mapping - new mapping data structure

173:     Level: developer

175:     Notes: Does not copy the indices, just keeps the pointer to the indices. The ISLocalToGlobalMappingDestroy()
176:     will free the space so it must be obtained with PetscMalloc() and it must not be freed elsewhere.

178:     Concepts: mapping^local to global

180: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate()
181: @*/
182: PetscErrorCode  ISLocalToGlobalMappingCreateNC(MPI_Comm cm,PetscInt n,const PetscInt indices[],ISLocalToGlobalMapping *mapping)
183: {

187:   if (n) {
189:   }
191:   *mapping = PETSC_NULL;
192: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
193:   VecInitializePackage(PETSC_NULL);
194: #endif
195:   if (IS_LTOGM_COOKIE == -1) {
196:     PetscLogClassRegister(&IS_LTOGM_COOKIE,"IS Local to global mapping");
197:   }

199:   PetscHeaderCreate(*mapping,_p_ISLocalToGlobalMapping,int,IS_LTOGM_COOKIE,0,"ISLocalToGlobalMapping",
200:                     cm,ISLocalToGlobalMappingDestroy,ISLocalToGlobalMappingView);
201:   PetscLogObjectMemory(*mapping,sizeof(struct _p_ISLocalToGlobalMapping)+n*sizeof(PetscInt));

203:   (*mapping)->n       = n;
204:   (*mapping)->indices = (PetscInt*)indices;

206:   /*
207:       Do not create the global to local mapping. This is only created if 
208:      ISGlobalToLocalMapping() is called 
209:   */
210:   (*mapping)->globals = 0;
211:   return(0);
212: }

216: /*@
217:     ISLocalToGlobalMappingBlock - Creates a blocked index version of an 
218:        ISLocalToGlobalMapping that is appropriate for MatSetLocalToGlobalMappingBlock()
219:        and VecSetLocalToGlobalMappingBlock().

221:     Not Collective, but communicator may have more than one process

223:     Input Parameters:
224: +    inmap - original point-wise mapping
225: -    bs - block size

227:     Output Parameter:
228: .   outmap - block based mapping; the indices are relative to BLOCKS, not individual vector or matrix entries.

230:     Level: advanced

232:     Concepts: mapping^local to global

234: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingCreateIS()
235: @*/
236: PetscErrorCode  ISLocalToGlobalMappingBlock(ISLocalToGlobalMapping inmap,PetscInt bs,ISLocalToGlobalMapping *outmap)
237: {
239:   PetscInt       *ii,i,n;

242:   if (bs > 1) {
243:     n    = inmap->n/bs;
244:     if (n*bs != inmap->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Pointwise mapping length is not divisible by block size");
245:     PetscMalloc(n*sizeof(PetscInt),&ii);
246:     for (i=0; i<n; i++) {
247:       ii[i] = inmap->indices[bs*i]/bs;
248:     }
249:     ISLocalToGlobalMappingCreate(inmap->comm,n,ii,outmap);
250:     PetscFree(ii);
251:   } else {
252:     *outmap = inmap;
253:     PetscObjectReference((PetscObject)inmap);
254:   }
255:   return(0);
256: }
257: 
260: /*@
261:    ISLocalToGlobalMappingDestroy - Destroys a mapping between a local (0 to n)
262:    ordering and a global parallel ordering.

264:    Note Collective

266:    Input Parameters:
267: .  mapping - mapping data structure

269:    Level: advanced

271: .seealso: ISLocalToGlobalMappingCreate()
272: @*/
273: PetscErrorCode  ISLocalToGlobalMappingDestroy(ISLocalToGlobalMapping mapping)
274: {
278:   if (--mapping->refct > 0) return(0);
279:   if (mapping->refct < 0) {
280:     SETERRQ(PETSC_ERR_PLIB,"Mapping already destroyed");
281:   }

283:   PetscFree(mapping->indices);
284:   PetscFree(mapping->globals);
285:   PetscHeaderDestroy(mapping);
286:   return(0);
287: }
288: 
291: /*@
292:     ISLocalToGlobalMappingApplyIS - Creates from an IS in the local numbering
293:     a new index set using the global numbering defined in an ISLocalToGlobalMapping
294:     context.

296:     Not collective

298:     Input Parameters:
299: +   mapping - mapping between local and global numbering
300: -   is - index set in local numbering

302:     Output Parameters:
303: .   newis - index set in global numbering

305:     Level: advanced

307:     Concepts: mapping^local to global

309: .seealso: ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingCreate(),
310:           ISLocalToGlobalMappingDestroy(), ISGlobalToLocalMappingApply()
311: @*/
312: PetscErrorCode  ISLocalToGlobalMappingApplyIS(ISLocalToGlobalMapping mapping,IS is,IS *newis)
313: {
315:   PetscInt            n,i,*idxin,*idxmap,*idxout,Nmax = mapping->n;


322:   ISGetLocalSize(is,&n);
323:   ISGetIndices(is,&idxin);
324:   idxmap = mapping->indices;
325: 
326:   PetscMalloc(n*sizeof(PetscInt),&idxout);
327:   for (i=0; i<n; i++) {
328:     if (idxin[i] >= Nmax) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"Local index %d too large %d (max) at %d",idxin[i],Nmax-1,i);
329:     idxout[i] = idxmap[idxin[i]];
330:   }
331:   ISRestoreIndices(is,&idxin);
332:   ISCreateGeneral(PETSC_COMM_SELF,n,idxout,newis);
333:   PetscFree(idxout);
334:   return(0);
335: }

337: /*MC
338:    ISLocalToGlobalMappingApply - Takes a list of integers in a local numbering
339:    and converts them to the global numbering.

341:    Not collective

343:    Input Parameters:
344: +  mapping - the local to global mapping context
345: .  N - number of integers
346: -  in - input indices in local numbering

348:    Output Parameter:
349: .  out - indices in global numbering

351:    Synopsis:
352:    PetscErrorCode ISLocalToGlobalMappingApply(ISLocalToGlobalMapping mapping,int N,int in[],int out[])

354:    Notes: 
355:    The in and out array parameters may be identical.

357:    Level: advanced

359: .seealso: ISLocalToGlobalMappingCreate(),ISLocalToGlobalMappingDestroy(), 
360:           ISLocalToGlobalMappingApplyIS(),AOCreateBasic(),AOApplicationToPetsc(),
361:           AOPetscToApplication(), ISGlobalToLocalMappingApply()

363:     Concepts: mapping^local to global

365: M*/

367: /* -----------------------------------------------------------------------------------------*/

371: /*
372:     Creates the global fields in the ISLocalToGlobalMapping structure
373: */
374: static PetscErrorCode ISGlobalToLocalMappingSetUp_Private(ISLocalToGlobalMapping mapping)
375: {
377:   PetscInt            i,*idx = mapping->indices,n = mapping->n,end,start,*globals;

380:   end   = 0;
381:   start = 100000000;

383:   for (i=0; i<n; i++) {
384:     if (idx[i] < 0) continue;
385:     if (idx[i] < start) start = idx[i];
386:     if (idx[i] > end)   end   = idx[i];
387:   }
388:   if (start > end) {start = 0; end = -1;}
389:   mapping->globalstart = start;
390:   mapping->globalend   = end;

392:   PetscMalloc((end-start+2)*sizeof(PetscInt),&globals);
393:   mapping->globals = globals;
394:   for (i=0; i<end-start+1; i++) {
395:     globals[i] = -1;
396:   }
397:   for (i=0; i<n; i++) {
398:     if (idx[i] < 0) continue;
399:     globals[idx[i] - start] = i;
400:   }

402:   PetscLogObjectMemory(mapping,(end-start+1)*sizeof(PetscInt));
403:   return(0);
404: }

408: /*@
409:     ISGlobalToLocalMappingApply - Provides the local numbering for a list of integers
410:     specified with a global numbering.

412:     Not collective

414:     Input Parameters:
415: +   mapping - mapping between local and global numbering
416: .   type - IS_GTOLM_MASK - replaces global indices with no local value with -1
417:            IS_GTOLM_DROP - drops the indices with no local value from the output list
418: .   n - number of global indices to map
419: -   idx - global indices to map

421:     Output Parameters:
422: +   nout - number of indices in output array (if type == IS_GTOLM_MASK then nout = n)
423: -   idxout - local index of each global index, one must pass in an array long enough 
424:              to hold all the indices. You can call ISGlobalToLocalMappingApply() with 
425:              idxout == PETSC_NULL to determine the required length (returned in nout)
426:              and then allocate the required space and call ISGlobalToLocalMappingApply()
427:              a second time to set the values.

429:     Notes:
430:     Either nout or idxout may be PETSC_NULL. idx and idxout may be identical.

432:     This is not scalable in memory usage. Each processor requires O(Nglobal) size 
433:     array to compute these.

435:     Level: advanced

437:     Concepts: mapping^global to local

439: .seealso: ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingCreate(),
440:           ISLocalToGlobalMappingDestroy()
441: @*/
442: PetscErrorCode  ISGlobalToLocalMappingApply(ISLocalToGlobalMapping mapping,ISGlobalToLocalMappingType type,
443:                                   PetscInt n,const PetscInt idx[],PetscInt *nout,PetscInt idxout[])
444: {
445:   PetscInt i,*globals,nf = 0,tmp,start,end;

449:   if (!mapping->globals) {
450:     ISGlobalToLocalMappingSetUp_Private(mapping);
451:   }
452:   globals = mapping->globals;
453:   start   = mapping->globalstart;
454:   end     = mapping->globalend;

456:   if (type == IS_GTOLM_MASK) {
457:     if (idxout) {
458:       for (i=0; i<n; i++) {
459:         if (idx[i] < 0) idxout[i] = idx[i];
460:         else if (idx[i] < start) idxout[i] = -1;
461:         else if (idx[i] > end)   idxout[i] = -1;
462:         else                     idxout[i] = globals[idx[i] - start];
463:       }
464:     }
465:     if (nout) *nout = n;
466:   } else {
467:     if (idxout) {
468:       for (i=0; i<n; i++) {
469:         if (idx[i] < 0) continue;
470:         if (idx[i] < start) continue;
471:         if (idx[i] > end) continue;
472:         tmp = globals[idx[i] - start];
473:         if (tmp < 0) continue;
474:         idxout[nf++] = tmp;
475:       }
476:     } else {
477:       for (i=0; i<n; i++) {
478:         if (idx[i] < 0) continue;
479:         if (idx[i] < start) continue;
480:         if (idx[i] > end) continue;
481:         tmp = globals[idx[i] - start];
482:         if (tmp < 0) continue;
483:         nf++;
484:       }
485:     }
486:     if (nout) *nout = nf;
487:   }

489:   return(0);
490: }

494: /*@C
495:     ISLocalToGlobalMappingGetInfo - Gets the neighbor information for each processor and 
496:      each index shared by more than one processor 

498:     Collective on ISLocalToGlobalMapping

500:     Input Parameters:
501: .   mapping - the mapping from local to global indexing

503:     Output Parameter:
504: +   nproc - number of processors that are connected to this one
505: .   proc - neighboring processors
506: .   numproc - number of indices for each subdomain (processor)
507: -   indices - indices of local nodes shared with neighbor (sorted by global numbering)

509:     Level: advanced

511:     Concepts: mapping^local to global

513:     Fortran Usage: 
514: $        ISLocalToGlobalMpngGetInfoSize(ISLocalToGlobalMapping,PetscInt nproc,PetscInt numprocmax,ierr) followed by 
515: $        ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping,PetscInt nproc, PetscInt procs[nproc],PetscInt numprocs[nproc],
516:           PetscInt indices[nproc][numprocmax],ierr)
517:         There is no ISLocalToGlobalMappingRestoreInfo() in Fortran. You must make sure that procs[], numprocs[] and 
518:         indices[][] are large enough arrays, either by allocating them dynamically or defining static ones large enough.


521: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
522:           ISLocalToGlobalMappingRestoreInfo()
523: @*/
524: PetscErrorCode  ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
525: {
527:   PetscMPIInt    size,rank,tag1,tag2,tag3,*len,*source,imdex;
528:   PetscInt       i,n = mapping->n,Ng,ng,max = 0,*lindices = mapping->indices;
529:   PetscInt       *nprocs,*owner,nsends,*sends,j,*starts,nmax,nrecvs,*recvs,proc;
530:   PetscInt       cnt,scale,*ownedsenders,*nownedsenders,rstart,nowned;
531:   PetscInt       node,nownedm,nt,*sends2,nsends2,*starts2,*lens2,*dest,nrecvs2,*starts3,*recvs2,k,*bprocs,*tmp;
532:   PetscInt       first_procs,first_numprocs,*first_indices;
533:   MPI_Request    *recv_waits,*send_waits;
534:   MPI_Status     recv_status,*send_status,*recv_statuses;
535:   MPI_Comm       comm = mapping->comm;
536:   PetscTruth     debug = PETSC_FALSE;

539:   MPI_Comm_size(comm,&size);
540:   MPI_Comm_rank(comm,&rank);
541:   if (size == 1) {
542:     *nproc         = 0;
543:     *procs         = PETSC_NULL;
544:     PetscMalloc(sizeof(PetscInt),numprocs);
545:     (*numprocs)[0] = 0;
546:     PetscMalloc(sizeof(PetscInt*),indices);
547:     (*indices)[0]  = PETSC_NULL;
548:     return(0);
549:   }

551:   PetscOptionsHasName(PETSC_NULL,"-islocaltoglobalmappinggetinfo_debug",&debug);

553:   /*
554:     Notes on ISLocalToGlobalMappingGetInfo

556:     globally owned node - the nodes that have been assigned to this processor in global
557:            numbering, just for this routine.

559:     nontrivial globally owned node - node assigned to this processor that is on a subdomain
560:            boundary (i.e. is has more than one local owner)

562:     locally owned node - node that exists on this processors subdomain

564:     nontrivial locally owned node - node that is not in the interior (i.e. has more than one
565:            local subdomain
566:   */
567:   PetscObjectGetNewTag((PetscObject)mapping,&tag1);
568:   PetscObjectGetNewTag((PetscObject)mapping,&tag2);
569:   PetscObjectGetNewTag((PetscObject)mapping,&tag3);

571:   for (i=0; i<n; i++) {
572:     if (lindices[i] > max) max = lindices[i];
573:   }
574:   MPI_Allreduce(&max,&Ng,1,MPIU_INT,MPI_MAX,comm);
575:   Ng++;
576:   MPI_Comm_size(comm,&size);
577:   MPI_Comm_rank(comm,&rank);
578:   scale  = Ng/size + 1;
579:   ng     = scale; if (rank == size-1) ng = Ng - scale*(size-1); ng = PetscMax(1,ng);
580:   rstart = scale*rank;

582:   /* determine ownership ranges of global indices */
583:   PetscMalloc(2*size*sizeof(PetscInt),&nprocs);
584:   PetscMemzero(nprocs,2*size*sizeof(PetscInt));

586:   /* determine owners of each local node  */
587:   PetscMalloc(n*sizeof(PetscInt),&owner);
588:   for (i=0; i<n; i++) {
589:     proc             = lindices[i]/scale; /* processor that globally owns this index */
590:     nprocs[2*proc+1] = 1;                 /* processor globally owns at least one of ours */
591:     owner[i]         = proc;
592:     nprocs[2*proc]++;                     /* count of how many that processor globally owns of ours */
593:   }
594:   nsends = 0; for (i=0; i<size; i++) nsends += nprocs[2*i+1];
595:   PetscInfo1(0,"Number of global owners for my local data %d\n",nsends);

597:   /* inform other processors of number of messages and max length*/
598:   PetscMaxSum(comm,nprocs,&nmax,&nrecvs);
599:   PetscInfo1(0,"Number of local owners for my global data %d\n",nrecvs);

601:   /* post receives for owned rows */
602:   PetscMalloc((2*nrecvs+1)*(nmax+1)*sizeof(PetscInt),&recvs);
603:   PetscMalloc((nrecvs+1)*sizeof(MPI_Request),&recv_waits);
604:   for (i=0; i<nrecvs; i++) {
605:     MPI_Irecv(recvs+2*nmax*i,2*nmax,MPIU_INT,MPI_ANY_SOURCE,tag1,comm,recv_waits+i);
606:   }

608:   /* pack messages containing lists of local nodes to owners */
609:   PetscMalloc((2*n+1)*sizeof(PetscInt),&sends);
610:   PetscMalloc((size+1)*sizeof(PetscInt),&starts);
611:   starts[0]  = 0;
612:   for (i=1; i<size; i++) { starts[i] = starts[i-1] + 2*nprocs[2*i-2];}
613:   for (i=0; i<n; i++) {
614:     sends[starts[owner[i]]++] = lindices[i];
615:     sends[starts[owner[i]]++] = i;
616:   }
617:   PetscFree(owner);
618:   starts[0]  = 0;
619:   for (i=1; i<size; i++) { starts[i] = starts[i-1] + 2*nprocs[2*i-2];}

621:   /* send the messages */
622:   PetscMalloc((nsends+1)*sizeof(MPI_Request),&send_waits);
623:   PetscMalloc((nsends+1)*sizeof(PetscInt),&dest);
624:   cnt = 0;
625:   for (i=0; i<size; i++) {
626:     if (nprocs[2*i]) {
627:       MPI_Isend(sends+starts[i],2*nprocs[2*i],MPIU_INT,i,tag1,comm,send_waits+cnt);
628:       dest[cnt] = i;
629:       cnt++;
630:     }
631:   }
632:   PetscFree(starts);

634:   /* wait on receives */
635:   PetscMalloc((nrecvs+1)*sizeof(PetscMPIInt),&source);
636:   PetscMalloc((nrecvs+1)*sizeof(PetscMPIInt),&len);
637:   cnt  = nrecvs;
638:   PetscMalloc((ng+1)*sizeof(PetscInt),&nownedsenders);
639:   PetscMemzero(nownedsenders,ng*sizeof(PetscInt));
640:   while (cnt) {
641:     MPI_Waitany(nrecvs,recv_waits,&imdex,&recv_status);
642:     /* unpack receives into our local space */
643:     MPI_Get_count(&recv_status,MPIU_INT,&len[imdex]);
644:     source[imdex]  = recv_status.MPI_SOURCE;
645:     len[imdex]     = len[imdex]/2;
646:     /* count how many local owners for each of my global owned indices */
647:     for (i=0; i<len[imdex]; i++) nownedsenders[recvs[2*imdex*nmax+2*i]-rstart]++;
648:     cnt--;
649:   }
650:   PetscFree(recv_waits);

652:   /* count how many globally owned indices are on an edge multiplied by how many processors own them. */
653:   nowned  = 0;
654:   nownedm = 0;
655:   for (i=0; i<ng; i++) {
656:     if (nownedsenders[i] > 1) {nownedm += nownedsenders[i]; nowned++;}
657:   }

659:   /* create single array to contain rank of all local owners of each globally owned index */
660:   PetscMalloc((nownedm+1)*sizeof(PetscInt),&ownedsenders);
661:   PetscMalloc((ng+1)*sizeof(PetscInt),&starts);
662:   starts[0] = 0;
663:   for (i=1; i<ng; i++) {
664:     if (nownedsenders[i-1] > 1) starts[i] = starts[i-1] + nownedsenders[i-1];
665:     else starts[i] = starts[i-1];
666:   }

668:   /* for each nontrival globally owned node list all arriving processors */
669:   for (i=0; i<nrecvs; i++) {
670:     for (j=0; j<len[i]; j++) {
671:       node = recvs[2*i*nmax+2*j]-rstart;
672:       if (nownedsenders[node] > 1) {
673:         ownedsenders[starts[node]++] = source[i];
674:       }
675:     }
676:   }

678:   if (debug) { /* -----------------------------------  */
679:     starts[0]    = 0;
680:     for (i=1; i<ng; i++) {
681:       if (nownedsenders[i-1] > 1) starts[i] = starts[i-1] + nownedsenders[i-1];
682:       else starts[i] = starts[i-1];
683:     }
684:     for (i=0; i<ng; i++) {
685:       if (nownedsenders[i] > 1) {
686:         PetscSynchronizedPrintf(comm,"[%d] global node %d local owner processors: ",rank,i+rstart);
687:         for (j=0; j<nownedsenders[i]; j++) {
688:           PetscSynchronizedPrintf(comm,"%d ",ownedsenders[starts[i]+j]);
689:         }
690:         PetscSynchronizedPrintf(comm,"\n");
691:       }
692:     }
693:     PetscSynchronizedFlush(comm);
694:   }/* -----------------------------------  */

696:   /* wait on original sends */
697:   if (nsends) {
698:     PetscMalloc(nsends*sizeof(MPI_Status),&send_status);
699:     MPI_Waitall(nsends,send_waits,send_status);
700:     PetscFree(send_status);
701:   }
702:   PetscFree(send_waits);
703:   PetscFree(sends);
704:   PetscFree(nprocs);

706:   /* pack messages to send back to local owners */
707:   starts[0]    = 0;
708:   for (i=1; i<ng; i++) {
709:     if (nownedsenders[i-1] > 1) starts[i] = starts[i-1] + nownedsenders[i-1];
710:     else starts[i] = starts[i-1];
711:   }
712:   nsends2 = nrecvs;
713:   PetscMalloc((nsends2+1)*sizeof(PetscInt),&nprocs); /* length of each message */
714:   for (i=0; i<nrecvs; i++) {
715:     nprocs[i] = 1;
716:     for (j=0; j<len[i]; j++) {
717:       node = recvs[2*i*nmax+2*j]-rstart;
718:       if (nownedsenders[node] > 1) {
719:         nprocs[i] += 2 + nownedsenders[node];
720:       }
721:     }
722:   }
723:   nt = 0; for (i=0; i<nsends2; i++) nt += nprocs[i];
724:   PetscMalloc((nt+1)*sizeof(PetscInt),&sends2);
725:   PetscMalloc((nsends2+1)*sizeof(PetscInt),&starts2);
726:   starts2[0] = 0; for (i=1; i<nsends2; i++) starts2[i] = starts2[i-1] + nprocs[i-1];
727:   /*
728:      Each message is 1 + nprocs[i] long, and consists of 
729:        (0) the number of nodes being sent back 
730:        (1) the local node number,
731:        (2) the number of processors sharing it,
732:        (3) the processors sharing it
733:   */
734:   for (i=0; i<nsends2; i++) {
735:     cnt = 1;
736:     sends2[starts2[i]] = 0;
737:     for (j=0; j<len[i]; j++) {
738:       node = recvs[2*i*nmax+2*j]-rstart;
739:       if (nownedsenders[node] > 1) {
740:         sends2[starts2[i]]++;
741:         sends2[starts2[i]+cnt++] = recvs[2*i*nmax+2*j+1];
742:         sends2[starts2[i]+cnt++] = nownedsenders[node];
743:         PetscMemcpy(&sends2[starts2[i]+cnt],&ownedsenders[starts[node]],nownedsenders[node]*sizeof(PetscInt));
744:         cnt += nownedsenders[node];
745:       }
746:     }
747:   }

749:   /* receive the message lengths */
750:   nrecvs2 = nsends;
751:   PetscMalloc((nrecvs2+1)*sizeof(PetscInt),&lens2);
752:   PetscMalloc((nrecvs2+1)*sizeof(PetscInt),&starts3);
753:   PetscMalloc((nrecvs2+1)*sizeof(MPI_Request),&recv_waits);
754:   for (i=0; i<nrecvs2; i++) {
755:     MPI_Irecv(&lens2[i],1,MPIU_INT,dest[i],tag2,comm,recv_waits+i);
756:   }

758:   /* send the message lengths */
759:   for (i=0; i<nsends2; i++) {
760:     MPI_Send(&nprocs[i],1,MPIU_INT,source[i],tag2,comm);
761:   }

763:   /* wait on receives of lens */
764:   if (nrecvs2) {
765:     PetscMalloc(nrecvs2*sizeof(MPI_Status),&recv_statuses);
766:     MPI_Waitall(nrecvs2,recv_waits,recv_statuses);
767:     PetscFree(recv_statuses);
768:   }
769:   PetscFree(recv_waits);

771:   starts3[0] = 0;
772:   nt         = 0;
773:   for (i=0; i<nrecvs2-1; i++) {
774:     starts3[i+1] = starts3[i] + lens2[i];
775:     nt          += lens2[i];
776:   }
777:   nt += lens2[nrecvs2-1];

779:   PetscMalloc((nt+1)*sizeof(PetscInt),&recvs2);
780:   PetscMalloc((nrecvs2+1)*sizeof(MPI_Request),&recv_waits);
781:   for (i=0; i<nrecvs2; i++) {
782:     MPI_Irecv(recvs2+starts3[i],lens2[i],MPIU_INT,dest[i],tag3,comm,recv_waits+i);
783:   }
784: 
785:   /* send the messages */
786:   PetscMalloc((nsends2+1)*sizeof(MPI_Request),&send_waits);
787:   for (i=0; i<nsends2; i++) {
788:     MPI_Isend(sends2+starts2[i],nprocs[i],MPIU_INT,source[i],tag3,comm,send_waits+i);
789:   }

791:   /* wait on receives */
792:   if (nrecvs2) {
793:     PetscMalloc(nrecvs2*sizeof(MPI_Status),&recv_statuses);
794:     MPI_Waitall(nrecvs2,recv_waits,recv_statuses);
795:     PetscFree(recv_statuses);
796:   }
797:   PetscFree(recv_waits);
798:   PetscFree(nprocs);

800:   if (debug) { /* -----------------------------------  */
801:     cnt = 0;
802:     for (i=0; i<nrecvs2; i++) {
803:       nt = recvs2[cnt++];
804:       for (j=0; j<nt; j++) {
805:         PetscSynchronizedPrintf(comm,"[%d] local node %d number of subdomains %d: ",rank,recvs2[cnt],recvs2[cnt+1]);
806:         for (k=0; k<recvs2[cnt+1]; k++) {
807:           PetscSynchronizedPrintf(comm,"%d ",recvs2[cnt+2+k]);
808:         }
809:         cnt += 2 + recvs2[cnt+1];
810:         PetscSynchronizedPrintf(comm,"\n");
811:       }
812:     }
813:     PetscSynchronizedFlush(comm);
814:   } /* -----------------------------------  */

816:   /* count number subdomains for each local node */
817:   PetscMalloc(size*sizeof(PetscInt),&nprocs);
818:   PetscMemzero(nprocs,size*sizeof(PetscInt));
819:   cnt  = 0;
820:   for (i=0; i<nrecvs2; i++) {
821:     nt = recvs2[cnt++];
822:     for (j=0; j<nt; j++) {
823:       for (k=0; k<recvs2[cnt+1]; k++) {
824:         nprocs[recvs2[cnt+2+k]]++;
825:       }
826:       cnt += 2 + recvs2[cnt+1];
827:     }
828:   }
829:   nt = 0; for (i=0; i<size; i++) nt += (nprocs[i] > 0);
830:   *nproc    = nt;
831:   PetscMalloc((nt+1)*sizeof(PetscInt),procs);
832:   PetscMalloc((nt+1)*sizeof(PetscInt),numprocs);
833:   PetscMalloc((nt+1)*sizeof(PetscInt*),indices);
834:   PetscMalloc(size*sizeof(PetscInt),&bprocs);
835:   cnt       = 0;
836:   for (i=0; i<size; i++) {
837:     if (nprocs[i] > 0) {
838:       bprocs[i]        = cnt;
839:       (*procs)[cnt]    = i;
840:       (*numprocs)[cnt] = nprocs[i];
841:       PetscMalloc(nprocs[i]*sizeof(PetscInt),&(*indices)[cnt]);
842:       cnt++;
843:     }
844:   }

846:   /* make the list of subdomains for each nontrivial local node */
847:   PetscMemzero(*numprocs,nt*sizeof(PetscInt));
848:   cnt  = 0;
849:   for (i=0; i<nrecvs2; i++) {
850:     nt = recvs2[cnt++];
851:     for (j=0; j<nt; j++) {
852:       for (k=0; k<recvs2[cnt+1]; k++) {
853:         (*indices)[bprocs[recvs2[cnt+2+k]]][(*numprocs)[bprocs[recvs2[cnt+2+k]]]++] = recvs2[cnt];
854:       }
855:       cnt += 2 + recvs2[cnt+1];
856:     }
857:   }
858:   PetscFree(bprocs);
859:   PetscFree(recvs2);

861:   /* sort the node indexing by their global numbers */
862:   nt = *nproc;
863:   for (i=0; i<nt; i++) {
864:     PetscMalloc(((*numprocs)[i])*sizeof(PetscInt),&tmp);
865:     for (j=0; j<(*numprocs)[i]; j++) {
866:       tmp[j] = lindices[(*indices)[i][j]];
867:     }
868:     PetscSortIntWithArray((*numprocs)[i],tmp,(*indices)[i]);
869:     PetscFree(tmp);
870:   }

872:   if (debug) { /* -----------------------------------  */
873:     nt = *nproc;
874:     for (i=0; i<nt; i++) {
875:       PetscSynchronizedPrintf(comm,"[%d] subdomain %d number of indices %d: ",rank,(*procs)[i],(*numprocs)[i]);
876:       for (j=0; j<(*numprocs)[i]; j++) {
877:         PetscSynchronizedPrintf(comm,"%d ",(*indices)[i][j]);
878:       }
879:       PetscSynchronizedPrintf(comm,"\n");
880:     }
881:     PetscSynchronizedFlush(comm);
882:   } /* -----------------------------------  */

884:   /* wait on sends */
885:   if (nsends2) {
886:     PetscMalloc(nsends2*sizeof(MPI_Status),&send_status);
887:     MPI_Waitall(nsends2,send_waits,send_status);
888:     PetscFree(send_status);
889:   }

891:   PetscFree(starts3);
892:   PetscFree(dest);
893:   PetscFree(send_waits);

895:   PetscFree(nownedsenders);
896:   PetscFree(ownedsenders);
897:   PetscFree(starts);
898:   PetscFree(starts2);
899:   PetscFree(lens2);

901:   PetscFree(source);
902:   PetscFree(len);
903:   PetscFree(recvs);
904:   PetscFree(nprocs);
905:   PetscFree(sends2);

907:   /* put the information about myself as the first entry in the list */
908:   first_procs    = (*procs)[0];
909:   first_numprocs = (*numprocs)[0];
910:   first_indices  = (*indices)[0];
911:   for (i=0; i<*nproc; i++) {
912:     if ((*procs)[i] == rank) {
913:       (*procs)[0]    = (*procs)[i];
914:       (*numprocs)[0] = (*numprocs)[i];
915:       (*indices)[0]  = (*indices)[i];
916:       (*procs)[i]    = first_procs;
917:       (*numprocs)[i] = first_numprocs;
918:       (*indices)[i]  = first_indices;
919:       break;
920:     }
921:   }
922:   return(0);
923: }

927: /*@C
928:     ISLocalToGlobalMappingRestoreInfo - Frees the memory allocated by ISLocalToGlobalMappingGetInfo()

930:     Collective on ISLocalToGlobalMapping

932:     Input Parameters:
933: .   mapping - the mapping from local to global indexing

935:     Output Parameter:
936: +   nproc - number of processors that are connected to this one
937: .   proc - neighboring processors
938: .   numproc - number of indices for each processor
939: -   indices - indices of local nodes shared with neighbor (sorted by global numbering)

941:     Level: advanced

943: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
944:           ISLocalToGlobalMappingGetInfo()
945: @*/
946: PetscErrorCode  ISLocalToGlobalMappingRestoreInfo(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
947: {
949:   PetscInt i;

952:   PetscFree(*procs);
953:   PetscFree(*numprocs);
954:   if (*indices) {
955:     PetscFree((*indices)[0]);
956:     for (i=1; i<*nproc; i++) {
957:       PetscFree((*indices)[i]);
958:     }
959:     PetscFree(*indices);
960:   }
961:   return(0);
962: }