Actual source code: isltog.c

petsc-3.11.2 2019-05-18
Report Typos and Errors

  2:  #include <petsc/private/isimpl.h>
  3:  #include <petsc/private/hashmapi.h>
  4:  #include <petscsf.h>
  5:  #include <petscviewer.h>

  7: PetscClassId IS_LTOGM_CLASSID;
  8: static PetscErrorCode  ISLocalToGlobalMappingGetBlockInfo_Private(ISLocalToGlobalMapping,PetscInt*,PetscInt**,PetscInt**,PetscInt***);

 10: typedef struct {
 11:   PetscInt *globals;
 12: } ISLocalToGlobalMapping_Basic;

 14: typedef struct {
 15:   PetscHMapI globalht;
 16: } ISLocalToGlobalMapping_Hash;


 19: PetscErrorCode ISGetPointRange(IS pointIS, PetscInt *pStart, PetscInt *pEnd, const PetscInt **points)
 20: {
 21:   PetscInt       numCells, step = 1;
 22:   PetscBool      isStride;

 26:   *pStart = 0;
 27:   *points = NULL;
 28:   ISGetLocalSize(pointIS, &numCells);
 29:   PetscObjectTypeCompare((PetscObject) pointIS, ISSTRIDE, &isStride);
 30:   if (isStride) {ISStrideGetInfo(pointIS, pStart, &step);}
 31:   *pEnd   = *pStart + numCells;
 32:   if (!isStride || step != 1) {ISGetIndices(pointIS, points);}
 33:   return(0);
 34: }

 36: PetscErrorCode ISRestorePointRange(IS pointIS, PetscInt *pStart, PetscInt *pEnd, const PetscInt **points)
 37: {
 38:   PetscInt       step = 1;
 39:   PetscBool      isStride;

 43:   PetscObjectTypeCompare((PetscObject) pointIS, ISSTRIDE, &isStride);
 44:   if (isStride) {ISStrideGetInfo(pointIS, pStart, &step);}
 45:   if (!isStride || step != 1) {ISGetIndices(pointIS, points);}
 46:   return(0);
 47: }

 49: PetscErrorCode ISGetPointSubrange(IS subpointIS, PetscInt pStart, PetscInt pEnd, const PetscInt *points)
 50: {

 54:   if (points) {
 55:     ISSetType(subpointIS, ISGENERAL);
 56:     ISGeneralSetIndices(subpointIS, pEnd-pStart, &points[pStart], PETSC_USE_POINTER);
 57:   } else {
 58:     ISSetType(subpointIS, ISSTRIDE);
 59:     ISStrideSetStride(subpointIS, pEnd-pStart, pStart, 1);
 60:   }
 61:   return(0);
 62: }

 64: /* -----------------------------------------------------------------------------------------*/

 66: /*
 67:     Creates the global mapping information in the ISLocalToGlobalMapping structure

 69:     If the user has not selected how to handle the global to local mapping then use HASH for "large" problems
 70: */
 71: static PetscErrorCode ISGlobalToLocalMappingSetUp(ISLocalToGlobalMapping mapping)
 72: {
 73:   PetscInt       i,*idx = mapping->indices,n = mapping->n,end,start;

 77:   if (mapping->data) return(0);
 78:   end   = 0;
 79:   start = PETSC_MAX_INT;

 81:   for (i=0; i<n; i++) {
 82:     if (idx[i] < 0) continue;
 83:     if (idx[i] < start) start = idx[i];
 84:     if (idx[i] > end)   end   = idx[i];
 85:   }
 86:   if (start > end) {start = 0; end = -1;}
 87:   mapping->globalstart = start;
 88:   mapping->globalend   = end;
 89:   if (!((PetscObject)mapping)->type_name) {
 90:     if ((end - start) > PetscMax(4*n,1000000)) {
 91:       ISLocalToGlobalMappingSetType(mapping,ISLOCALTOGLOBALMAPPINGHASH);
 92:     } else {
 93:       ISLocalToGlobalMappingSetType(mapping,ISLOCALTOGLOBALMAPPINGBASIC);
 94:     }
 95:   }
 96:   (*mapping->ops->globaltolocalmappingsetup)(mapping);
 97:   return(0);
 98: }

100: static PetscErrorCode ISGlobalToLocalMappingSetUp_Basic(ISLocalToGlobalMapping mapping)
101: {
102:   PetscErrorCode              ierr;
103:   PetscInt                    i,*idx = mapping->indices,n = mapping->n,end,start,*globals;
104:   ISLocalToGlobalMapping_Basic *map;

107:   start            = mapping->globalstart;
108:   end              = mapping->globalend;
109:   PetscNew(&map);
110:   PetscMalloc1(end-start+2,&globals);
111:   map->globals     = globals;
112:   for (i=0; i<end-start+1; i++) globals[i] = -1;
113:   for (i=0; i<n; i++) {
114:     if (idx[i] < 0) continue;
115:     globals[idx[i] - start] = i;
116:   }
117:   mapping->data = (void*)map;
118:   PetscLogObjectMemory((PetscObject)mapping,(end-start+1)*sizeof(PetscInt));
119:   return(0);
120: }

122: static PetscErrorCode ISGlobalToLocalMappingSetUp_Hash(ISLocalToGlobalMapping mapping)
123: {
124:   PetscErrorCode              ierr;
125:   PetscInt                    i,*idx = mapping->indices,n = mapping->n;
126:   ISLocalToGlobalMapping_Hash *map;

129:   PetscNew(&map);
130:   PetscHMapICreate(&map->globalht);
131:   for (i=0; i<n; i++ ) {
132:     if (idx[i] < 0) continue;
133:     PetscHMapISet(map->globalht,idx[i],i);
134:   }
135:   mapping->data = (void*)map;
136:   PetscLogObjectMemory((PetscObject)mapping,2*n*sizeof(PetscInt));
137:   return(0);
138: }

140: static PetscErrorCode ISLocalToGlobalMappingDestroy_Basic(ISLocalToGlobalMapping mapping)
141: {
142:   PetscErrorCode              ierr;
143:   ISLocalToGlobalMapping_Basic *map  = (ISLocalToGlobalMapping_Basic *)mapping->data;

146:   if (!map) return(0);
147:   PetscFree(map->globals);
148:   PetscFree(mapping->data);
149:   return(0);
150: }

152: static PetscErrorCode ISLocalToGlobalMappingDestroy_Hash(ISLocalToGlobalMapping mapping)
153: {
154:   PetscErrorCode              ierr;
155:   ISLocalToGlobalMapping_Hash *map  = (ISLocalToGlobalMapping_Hash*)mapping->data;

158:   if (!map) return(0);
159:   PetscHMapIDestroy(&map->globalht);
160:   PetscFree(mapping->data);
161:   return(0);
162: }

164: #define GTOLTYPE _Basic
165: #define GTOLNAME _Basic
166: #define GTOLBS mapping->bs
167: #define GTOL(g, local) do {                  \
168:     local = map->globals[g/bs - start];      \
169:     local = bs*local + (g % bs);             \
170:   } while (0)

172:  #include <../src/vec/is/utils/isltog.h>

174: #define GTOLTYPE _Basic
175: #define GTOLNAME Block_Basic
176: #define GTOLBS 1
177: #define GTOL(g, local) do {                  \
178:     local = map->globals[g - start];         \
179:   } while (0)
180:  #include <../src/vec/is/utils/isltog.h>

182: #define GTOLTYPE _Hash
183: #define GTOLNAME _Hash
184: #define GTOLBS mapping->bs
185: #define GTOL(g, local) do {                         \
186:     (void)PetscHMapIGet(map->globalht,g/bs,&local); \
187:     local = bs*local + (g % bs);                    \
188:    } while (0)
189:  #include <../src/vec/is/utils/isltog.h>

191: #define GTOLTYPE _Hash
192: #define GTOLNAME Block_Hash
193: #define GTOLBS 1
194: #define GTOL(g, local) do {                         \
195:     (void)PetscHMapIGet(map->globalht,g,&local);    \
196:   } while (0)
197:  #include <../src/vec/is/utils/isltog.h>

199: /*@
200:     ISLocalToGlobalMappingDuplicate - Duplicates the local to global mapping object

202:     Not Collective

204:     Input Parameter:
205: .   ltog - local to global mapping

207:     Output Parameter:
208: .   nltog - the duplicated local to global mapping

210:     Level: advanced

212:     Concepts: mapping^local to global

214: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate()
215: @*/
216: PetscErrorCode  ISLocalToGlobalMappingDuplicate(ISLocalToGlobalMapping ltog,ISLocalToGlobalMapping* nltog)
217: {

222:   ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)ltog),ltog->bs,ltog->n,ltog->indices,PETSC_COPY_VALUES,nltog);
223:   return(0);
224: }

226: /*@
227:     ISLocalToGlobalMappingGetSize - Gets the local size of a local to global mapping

229:     Not Collective

231:     Input Parameter:
232: .   ltog - local to global mapping

234:     Output Parameter:
235: .   n - the number of entries in the local mapping, ISLocalToGlobalMappingGetIndices() returns an array of this length

237:     Level: advanced

239:     Concepts: mapping^local to global

241: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate()
242: @*/
243: PetscErrorCode  ISLocalToGlobalMappingGetSize(ISLocalToGlobalMapping mapping,PetscInt *n)
244: {
248:   *n = mapping->bs*mapping->n;
249:   return(0);
250: }

252: /*@C
253:     ISLocalToGlobalMappingView - View a local to global mapping

255:     Not Collective

257:     Input Parameters:
258: +   ltog - local to global mapping
259: -   viewer - viewer

261:     Level: advanced

263:     Concepts: mapping^local to global

265: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate()
266: @*/
267: PetscErrorCode  ISLocalToGlobalMappingView(ISLocalToGlobalMapping mapping,PetscViewer viewer)
268: {
269:   PetscInt       i;
270:   PetscMPIInt    rank;
271:   PetscBool      iascii;

276:   if (!viewer) {
277:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mapping),&viewer);
278:   }

281:   MPI_Comm_rank(PetscObjectComm((PetscObject)mapping),&rank);
282:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
283:   if (iascii) {
284:     PetscObjectPrintClassNamePrefixType((PetscObject)mapping,viewer);
285:     PetscViewerASCIIPushSynchronized(viewer);
286:     for (i=0; i<mapping->n; i++) {
287:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,mapping->indices[i]);
288:     }
289:     PetscViewerFlush(viewer);
290:     PetscViewerASCIIPopSynchronized(viewer);
291:   }
292:   return(0);
293: }

295: /*@
296:     ISLocalToGlobalMappingCreateIS - Creates a mapping between a local (0 to n)
297:     ordering and a global parallel ordering.

299:     Not collective

301:     Input Parameter:
302: .   is - index set containing the global numbers for each local number

304:     Output Parameter:
305: .   mapping - new mapping data structure

307:     Notes:
308:     the block size of the IS determines the block size of the mapping
309:     Level: advanced

311:     Concepts: mapping^local to global

313: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingSetFromOptions()
314: @*/
315: PetscErrorCode  ISLocalToGlobalMappingCreateIS(IS is,ISLocalToGlobalMapping *mapping)
316: {
318:   PetscInt       n,bs;
319:   const PetscInt *indices;
320:   MPI_Comm       comm;
321:   PetscBool      isblock;


327:   PetscObjectGetComm((PetscObject)is,&comm);
328:   ISGetLocalSize(is,&n);
329:   PetscObjectTypeCompare((PetscObject)is,ISBLOCK,&isblock);
330:   if (!isblock) {
331:     ISGetIndices(is,&indices);
332:     ISLocalToGlobalMappingCreate(comm,1,n,indices,PETSC_COPY_VALUES,mapping);
333:     ISRestoreIndices(is,&indices);
334:   } else {
335:     ISGetBlockSize(is,&bs);
336:     ISBlockGetIndices(is,&indices);
337:     ISLocalToGlobalMappingCreate(comm,bs,n/bs,indices,PETSC_COPY_VALUES,mapping);
338:     ISBlockRestoreIndices(is,&indices);
339:   }
340:   return(0);
341: }

343: /*@C
344:     ISLocalToGlobalMappingCreateSF - Creates a mapping between a local (0 to n)
345:     ordering and a global parallel ordering.

347:     Collective

349:     Input Parameter:
350: +   sf - star forest mapping contiguous local indices to (rank, offset)
351: -   start - first global index on this process

353:     Output Parameter:
354: .   mapping - new mapping data structure

356:     Level: advanced

358:     Concepts: mapping^local to global

360: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingSetFromOptions()
361: @*/
362: PetscErrorCode ISLocalToGlobalMappingCreateSF(PetscSF sf,PetscInt start,ISLocalToGlobalMapping *mapping)
363: {
365:   PetscInt       i,maxlocal,nroots,nleaves,*globals,*ltog;
366:   const PetscInt *ilocal;
367:   MPI_Comm       comm;


373:   PetscObjectGetComm((PetscObject)sf,&comm);
374:   PetscSFGetGraph(sf,&nroots,&nleaves,&ilocal,NULL);
375:   if (ilocal) {
376:     for (i=0,maxlocal=0; i<nleaves; i++) maxlocal = PetscMax(maxlocal,ilocal[i]+1);
377:   }
378:   else maxlocal = nleaves;
379:   PetscMalloc1(nroots,&globals);
380:   PetscMalloc1(maxlocal,&ltog);
381:   for (i=0; i<nroots; i++) globals[i] = start + i;
382:   for (i=0; i<maxlocal; i++) ltog[i] = -1;
383:   PetscSFBcastBegin(sf,MPIU_INT,globals,ltog);
384:   PetscSFBcastEnd(sf,MPIU_INT,globals,ltog);
385:   ISLocalToGlobalMappingCreate(comm,1,maxlocal,ltog,PETSC_OWN_POINTER,mapping);
386:   PetscFree(globals);
387:   return(0);
388: }

390: /*@
391:     ISLocalToGlobalMappingSetBlockSize - Sets the blocksize of the mapping

393:     Not collective

395:     Input Parameters:
396: .   mapping - mapping data structure
397: .   bs - the blocksize

399:     Level: advanced

401:     Concepts: mapping^local to global

403: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS()
404: @*/
405: PetscErrorCode  ISLocalToGlobalMappingSetBlockSize(ISLocalToGlobalMapping mapping,PetscInt bs)
406: {
407:   PetscInt       *nid;
408:   const PetscInt *oid;
409:   PetscInt       i,cn,on,obs,nn;

414:   if (bs < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid block size %D",bs);
415:   if (bs == mapping->bs) return(0);
416:   on  = mapping->n;
417:   obs = mapping->bs;
418:   oid = mapping->indices;
419:   nn  = (on*obs)/bs;
420:   if ((on*obs)%bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Block size %D is inconsistent with block size %D and number of block indices %D",bs,obs,on);

422:   PetscMalloc1(nn,&nid);
423:   ISLocalToGlobalMappingGetIndices(mapping,&oid);
424:   for (i=0;i<nn;i++) {
425:     PetscInt j;
426:     for (j=0,cn=0;j<bs-1;j++) {
427:       if (oid[i*bs+j] < 0) { cn++; continue; }
428:       if (oid[i*bs+j] != oid[i*bs+j+1]-1) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Block sizes %D and %D are incompatible with the block indices: non consecutive indices %D %D",bs,obs,oid[i*bs+j],oid[i*bs+j+1]);
429:     }
430:     if (oid[i*bs+j] < 0) cn++;
431:     if (cn) {
432:       if (cn != bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Block sizes %D and %D are incompatible with the block indices: invalid number of negative entries in block %D",bs,obs,cn);
433:       nid[i] = -1;
434:     } else {
435:       nid[i] = oid[i*bs]/bs;
436:     }
437:   }
438:   ISLocalToGlobalMappingRestoreIndices(mapping,&oid);

440:   mapping->n           = nn;
441:   mapping->bs          = bs;
442:   PetscFree(mapping->indices);
443:   mapping->indices     = nid;
444:   mapping->globalstart = 0;
445:   mapping->globalend   = 0;

447:   /* reset the cached information */
448:   PetscFree(mapping->info_procs);
449:   PetscFree(mapping->info_numprocs);
450:   if (mapping->info_indices) {
451:     PetscInt i;

453:     PetscFree((mapping->info_indices)[0]);
454:     for (i=1; i<mapping->info_nproc; i++) {
455:       PetscFree(mapping->info_indices[i]);
456:     }
457:     PetscFree(mapping->info_indices);
458:   }
459:   mapping->info_cached = PETSC_FALSE;

461:   if (mapping->ops->destroy) {
462:     (*mapping->ops->destroy)(mapping);
463:   }
464:   return(0);
465: }

467: /*@
468:     ISLocalToGlobalMappingGetBlockSize - Gets the blocksize of the mapping
469:     ordering and a global parallel ordering.

471:     Not Collective

473:     Input Parameters:
474: .   mapping - mapping data structure

476:     Output Parameter:
477: .   bs - the blocksize

479:     Level: advanced

481:     Concepts: mapping^local to global

483: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS()
484: @*/
485: PetscErrorCode  ISLocalToGlobalMappingGetBlockSize(ISLocalToGlobalMapping mapping,PetscInt *bs)
486: {
489:   *bs = mapping->bs;
490:   return(0);
491: }

493: /*@
494:     ISLocalToGlobalMappingCreate - Creates a mapping between a local (0 to n)
495:     ordering and a global parallel ordering.

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

499:     Input Parameters:
500: +   comm - MPI communicator
501: .   bs - the block size
502: .   n - the number of local elements divided by the block size, or equivalently the number of block indices
503: .   indices - the global index for each local element, these do not need to be in increasing order (sorted), these values should not be scaled (i.e. multiplied) by the blocksize bs
504: -   mode - see PetscCopyMode

506:     Output Parameter:
507: .   mapping - new mapping data structure

509:     Notes:
510:     There is one integer value in indices per block and it represents the actual indices bs*idx + j, where j=0,..,bs-1

512:     For "small" problems when using ISGlobalToLocalMappingApply() and ISGlobalToLocalMappingApplyBlock(), the ISLocalToGlobalMappingType of ISLOCALTOGLOBALMAPPINGBASIC will be used;
513:     this uses more memory but is faster; this approach is not scalable for extremely large mappings. For large problems ISLOCALTOGLOBALMAPPINGHASH is used, this is scalable.
514:     Use ISLocalToGlobalMappingSetType() or call ISLocalToGlobalMappingSetFromOptions() with the option -islocaltoglobalmapping_type <basic,hash> to control which is used.

516:     Level: advanced

518:     Concepts: mapping^local to global

520: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingSetFromOptions(), ISLOCALTOGLOBALMAPPINGBASIC, ISLOCALTOGLOBALMAPPINGHASH
521:           ISLocalToGlobalMappingSetType(), ISLocalToGlobalMappingType
522: @*/
523: PetscErrorCode  ISLocalToGlobalMappingCreate(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscInt indices[],PetscCopyMode mode,ISLocalToGlobalMapping *mapping)
524: {
526:   PetscInt       *in;


532:   *mapping = NULL;
533:   ISInitializePackage();

535:   PetscHeaderCreate(*mapping,IS_LTOGM_CLASSID,"ISLocalToGlobalMapping","Local to global mapping","IS",
536:                            comm,ISLocalToGlobalMappingDestroy,ISLocalToGlobalMappingView);
537:   (*mapping)->n             = n;
538:   (*mapping)->bs            = bs;
539:   (*mapping)->info_cached   = PETSC_FALSE;
540:   (*mapping)->info_free     = PETSC_FALSE;
541:   (*mapping)->info_procs    = NULL;
542:   (*mapping)->info_numprocs = NULL;
543:   (*mapping)->info_indices  = NULL;
544:   (*mapping)->info_nodec    = NULL;
545:   (*mapping)->info_nodei    = NULL;

547:   (*mapping)->ops->globaltolocalmappingapply      = NULL;
548:   (*mapping)->ops->globaltolocalmappingapplyblock = NULL;
549:   (*mapping)->ops->destroy                        = NULL;
550:   if (mode == PETSC_COPY_VALUES) {
551:     PetscMalloc1(n,&in);
552:     PetscMemcpy(in,indices,n*sizeof(PetscInt));
553:     (*mapping)->indices = in;
554:     PetscLogObjectMemory((PetscObject)*mapping,n*sizeof(PetscInt));
555:   } else if (mode == PETSC_OWN_POINTER) {
556:     (*mapping)->indices = (PetscInt*)indices;
557:     PetscLogObjectMemory((PetscObject)*mapping,n*sizeof(PetscInt));
558:   }
559:   else SETERRQ(comm,PETSC_ERR_SUP,"Cannot currently use PETSC_USE_POINTER");
560:   return(0);
561: }

563: PetscFunctionList ISLocalToGlobalMappingList = NULL;


566: /*@
567:    ISLocalToGlobalMappingSetFromOptions - Set mapping options from the options database.

569:    Not collective

571:    Input Parameters:
572: .  mapping - mapping data structure

574:    Level: advanced

576: @*/
577: PetscErrorCode ISLocalToGlobalMappingSetFromOptions(ISLocalToGlobalMapping mapping)
578: {
579:   PetscErrorCode             ierr;
580:   char                       type[256];
581:   ISLocalToGlobalMappingType defaulttype = "Not set";
582:   PetscBool                  flg;

586:   ISLocalToGlobalMappingRegisterAll();
587:   PetscObjectOptionsBegin((PetscObject)mapping);
588:   PetscOptionsFList("-islocaltoglobalmapping_type","ISLocalToGlobalMapping method","ISLocalToGlobalMappingSetType",ISLocalToGlobalMappingList,(char*)(((PetscObject)mapping)->type_name) ? ((PetscObject)mapping)->type_name : defaulttype,type,256,&flg);
589:   if (flg) {
590:     ISLocalToGlobalMappingSetType(mapping,type);
591:   }
592:   PetscOptionsEnd();
593:   return(0);
594: }

596: /*@
597:    ISLocalToGlobalMappingDestroy - Destroys a mapping between a local (0 to n)
598:    ordering and a global parallel ordering.

600:    Note Collective

602:    Input Parameters:
603: .  mapping - mapping data structure

605:    Level: advanced

607: .seealso: ISLocalToGlobalMappingCreate()
608: @*/
609: PetscErrorCode  ISLocalToGlobalMappingDestroy(ISLocalToGlobalMapping *mapping)
610: {

614:   if (!*mapping) return(0);
616:   if (--((PetscObject)(*mapping))->refct > 0) {*mapping = 0;return(0);}
617:   PetscFree((*mapping)->indices);
618:   PetscFree((*mapping)->info_procs);
619:   PetscFree((*mapping)->info_numprocs);
620:   if ((*mapping)->info_indices) {
621:     PetscInt i;

623:     PetscFree(((*mapping)->info_indices)[0]);
624:     for (i=1; i<(*mapping)->info_nproc; i++) {
625:       PetscFree(((*mapping)->info_indices)[i]);
626:     }
627:     PetscFree((*mapping)->info_indices);
628:   }
629:   if ((*mapping)->info_nodei) {
630:     PetscFree(((*mapping)->info_nodei)[0]);
631:   }
632:   PetscFree((*mapping)->info_nodei);
633:   PetscFree((*mapping)->info_nodec);
634:   if ((*mapping)->ops->destroy) {
635:     (*(*mapping)->ops->destroy)(*mapping);
636:   }
637:   PetscHeaderDestroy(mapping);
638:   *mapping = 0;
639:   return(0);
640: }

642: /*@
643:     ISLocalToGlobalMappingApplyIS - Creates from an IS in the local numbering
644:     a new index set using the global numbering defined in an ISLocalToGlobalMapping
645:     context.

647:     Collective on is

649:     Input Parameters:
650: +   mapping - mapping between local and global numbering
651: -   is - index set in local numbering

653:     Output Parameters:
654: .   newis - index set in global numbering

656:     Notes:
657:     The output IS will have the same communicator of the input IS.

659:     Level: advanced

661:     Concepts: mapping^local to global

663: .seealso: ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingCreate(),
664:           ISLocalToGlobalMappingDestroy(), ISGlobalToLocalMappingApply()
665: @*/
666: PetscErrorCode  ISLocalToGlobalMappingApplyIS(ISLocalToGlobalMapping mapping,IS is,IS *newis)
667: {
669:   PetscInt       n,*idxout;
670:   const PetscInt *idxin;


677:   ISGetLocalSize(is,&n);
678:   ISGetIndices(is,&idxin);
679:   PetscMalloc1(n,&idxout);
680:   ISLocalToGlobalMappingApply(mapping,n,idxin,idxout);
681:   ISRestoreIndices(is,&idxin);
682:   ISCreateGeneral(PetscObjectComm((PetscObject)is),n,idxout,PETSC_OWN_POINTER,newis);
683:   return(0);
684: }

686: /*@
687:    ISLocalToGlobalMappingApply - Takes a list of integers in a local numbering
688:    and converts them to the global numbering.

690:    Not collective

692:    Input Parameters:
693: +  mapping - the local to global mapping context
694: .  N - number of integers
695: -  in - input indices in local numbering

697:    Output Parameter:
698: .  out - indices in global numbering

700:    Notes:
701:    The in and out array parameters may be identical.

703:    Level: advanced

705: .seealso: ISLocalToGlobalMappingApplyBlock(), ISLocalToGlobalMappingCreate(),ISLocalToGlobalMappingDestroy(),
706:           ISLocalToGlobalMappingApplyIS(),AOCreateBasic(),AOApplicationToPetsc(),
707:           AOPetscToApplication(), ISGlobalToLocalMappingApply()

709:     Concepts: mapping^local to global
710: @*/
711: PetscErrorCode ISLocalToGlobalMappingApply(ISLocalToGlobalMapping mapping,PetscInt N,const PetscInt in[],PetscInt out[])
712: {
713:   PetscInt i,bs,Nmax;

717:   bs   = mapping->bs;
718:   Nmax = bs*mapping->n;
719:   if (bs == 1) {
720:     const PetscInt *idx = mapping->indices;
721:     for (i=0; i<N; i++) {
722:       if (in[i] < 0) {
723:         out[i] = in[i];
724:         continue;
725:       }
726:       if (in[i] >= Nmax) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local index %D too large %D (max) at %D",in[i],Nmax-1,i);
727:       out[i] = idx[in[i]];
728:     }
729:   } else {
730:     const PetscInt *idx = mapping->indices;
731:     for (i=0; i<N; i++) {
732:       if (in[i] < 0) {
733:         out[i] = in[i];
734:         continue;
735:       }
736:       if (in[i] >= Nmax) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local index %D too large %D (max) at %D",in[i],Nmax-1,i);
737:       out[i] = idx[in[i]/bs]*bs + (in[i] % bs);
738:     }
739:   }
740:   return(0);
741: }

743: /*@
744:    ISLocalToGlobalMappingApplyBlock - Takes a list of integers in a local block numbering and converts them to the global block numbering

746:    Not collective

748:    Input Parameters:
749: +  mapping - the local to global mapping context
750: .  N - number of integers
751: -  in - input indices in local block numbering

753:    Output Parameter:
754: .  out - indices in global block numbering

756:    Notes:
757:    The in and out array parameters may be identical.

759:    Example:
760:      If the index values are {0,1,6,7} set with a call to ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,2,2,{0,3}) then the mapping applied to 0
761:      (the first block) would produce 0 and the mapping applied to 1 (the second block) would produce 3.

763:    Level: advanced

765: .seealso: ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingCreate(),ISLocalToGlobalMappingDestroy(),
766:           ISLocalToGlobalMappingApplyIS(),AOCreateBasic(),AOApplicationToPetsc(),
767:           AOPetscToApplication(), ISGlobalToLocalMappingApply()

769:     Concepts: mapping^local to global
770: @*/
771: PetscErrorCode ISLocalToGlobalMappingApplyBlock(ISLocalToGlobalMapping mapping,PetscInt N,const PetscInt in[],PetscInt out[])
772: {

776:   {
777:     PetscInt i,Nmax = mapping->n;
778:     const PetscInt *idx = mapping->indices;

780:     for (i=0; i<N; i++) {
781:       if (in[i] < 0) {
782:         out[i] = in[i];
783:         continue;
784:       }
785:       if (in[i] >= Nmax) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local block index %D too large %D (max) at %D",in[i],Nmax-1,i);
786:       out[i] = idx[in[i]];
787:     }
788:   }
789:   return(0);
790: }

792: /*@
793:     ISGlobalToLocalMappingApply - Provides the local numbering for a list of integers
794:     specified with a global numbering.

796:     Not collective

798:     Input Parameters:
799: +   mapping - mapping between local and global numbering
800: .   type - IS_GTOLM_MASK - replaces global indices with no local value with -1
801:            IS_GTOLM_DROP - drops the indices with no local value from the output list
802: .   n - number of global indices to map
803: -   idx - global indices to map

805:     Output Parameters:
806: +   nout - number of indices in output array (if type == IS_GTOLM_MASK then nout = n)
807: -   idxout - local index of each global index, one must pass in an array long enough
808:              to hold all the indices. You can call ISGlobalToLocalMappingApply() with
809:              idxout == NULL to determine the required length (returned in nout)
810:              and then allocate the required space and call ISGlobalToLocalMappingApply()
811:              a second time to set the values.

813:     Notes:
814:     Either nout or idxout may be NULL. idx and idxout may be identical.

816:     For "small" problems when using ISGlobalToLocalMappingApply() and ISGlobalToLocalMappingApplyBlock(), the ISLocalToGlobalMappingType of ISLOCALTOGLOBALMAPPINGBASIC will be used;
817:     this uses more memory but is faster; this approach is not scalable for extremely large mappings. For large problems ISLOCALTOGLOBALMAPPINGHASH is used, this is scalable.
818:     Use ISLocalToGlobalMappingSetType() or call ISLocalToGlobalMappingSetFromOptions() with the option -islocaltoglobalmapping_type <basic,hash> to control which is used.

820:     Level: advanced

822:     Developer Note: The manual page states that idx and idxout may be identical but the calling
823:        sequence declares idx as const so it cannot be the same as idxout.

825:     Concepts: mapping^global to local

827: .seealso: ISLocalToGlobalMappingApply(), ISGlobalToLocalMappingApplyBlock(), ISLocalToGlobalMappingCreate(),
828:           ISLocalToGlobalMappingDestroy()
829: @*/
830: PetscErrorCode  ISGlobalToLocalMappingApply(ISLocalToGlobalMapping mapping,ISGlobalToLocalMappingMode type,PetscInt n,const PetscInt idx[],PetscInt *nout,PetscInt idxout[])
831: {

836:   if (!mapping->data) {
837:     ISGlobalToLocalMappingSetUp(mapping);
838:   }
839:   (*mapping->ops->globaltolocalmappingapply)(mapping,type,n,idx,nout,idxout);
840:   return(0);
841: }

843: /*@
844:     ISGlobalToLocalMappingApplyIS - Creates from an IS in the global numbering
845:     a new index set using the local numbering defined in an ISLocalToGlobalMapping
846:     context.

848:     Not collective

850:     Input Parameters:
851: +   mapping - mapping between local and global numbering
852: .   type - IS_GTOLM_MASK - replaces global indices with no local value with -1
853:            IS_GTOLM_DROP - drops the indices with no local value from the output list
854: -   is - index set in global numbering

856:     Output Parameters:
857: .   newis - index set in local numbering

859:     Notes:
860:     The output IS will be sequential, as it encodes a purely local operation

862:     Level: advanced

864:     Concepts: mapping^local to global

866: .seealso: ISGlobalToLocalMappingApply(), ISLocalToGlobalMappingCreate(),
867:           ISLocalToGlobalMappingDestroy()
868: @*/
869: PetscErrorCode  ISGlobalToLocalMappingApplyIS(ISLocalToGlobalMapping mapping,ISGlobalToLocalMappingMode type,IS is,IS *newis)
870: {
872:   PetscInt       n,nout,*idxout;
873:   const PetscInt *idxin;


880:   ISGetLocalSize(is,&n);
881:   ISGetIndices(is,&idxin);
882:   if (type == IS_GTOLM_MASK) {
883:     PetscMalloc1(n,&idxout);
884:   } else {
885:     ISGlobalToLocalMappingApply(mapping,type,n,idxin,&nout,NULL);
886:     PetscMalloc1(nout,&idxout);
887:   }
888:   ISGlobalToLocalMappingApply(mapping,type,n,idxin,&nout,idxout);
889:   ISRestoreIndices(is,&idxin);
890:   ISCreateGeneral(PETSC_COMM_SELF,nout,idxout,PETSC_OWN_POINTER,newis);
891:   return(0);
892: }

894: /*@
895:     ISGlobalToLocalMappingApplyBlock - Provides the local block numbering for a list of integers
896:     specified with a block global numbering.

898:     Not collective

900:     Input Parameters:
901: +   mapping - mapping between local and global numbering
902: .   type - IS_GTOLM_MASK - replaces global indices with no local value with -1
903:            IS_GTOLM_DROP - drops the indices with no local value from the output list
904: .   n - number of global indices to map
905: -   idx - global indices to map

907:     Output Parameters:
908: +   nout - number of indices in output array (if type == IS_GTOLM_MASK then nout = n)
909: -   idxout - local index of each global index, one must pass in an array long enough
910:              to hold all the indices. You can call ISGlobalToLocalMappingApplyBlock() with
911:              idxout == NULL to determine the required length (returned in nout)
912:              and then allocate the required space and call ISGlobalToLocalMappingApplyBlock()
913:              a second time to set the values.

915:     Notes:
916:     Either nout or idxout may be NULL. idx and idxout may be identical.

918:     For "small" problems when using ISGlobalToLocalMappingApply() and ISGlobalToLocalMappingApplyBlock(), the ISLocalToGlobalMappingType of ISLOCALTOGLOBALMAPPINGBASIC will be used;
919:     this uses more memory but is faster; this approach is not scalable for extremely large mappings. For large problems ISLOCALTOGLOBALMAPPINGHASH is used, this is scalable.
920:     Use ISLocalToGlobalMappingSetType() or call ISLocalToGlobalMappingSetFromOptions() with the option -islocaltoglobalmapping_type <basic,hash> to control which is used.

922:     Level: advanced

924:     Developer Note: The manual page states that idx and idxout may be identical but the calling
925:        sequence declares idx as const so it cannot be the same as idxout.

927:     Concepts: mapping^global to local

929: .seealso: ISLocalToGlobalMappingApply(), ISGlobalToLocalMappingApply(), ISLocalToGlobalMappingCreate(),
930:           ISLocalToGlobalMappingDestroy()
931: @*/
932: PetscErrorCode  ISGlobalToLocalMappingApplyBlock(ISLocalToGlobalMapping mapping,ISGlobalToLocalMappingMode type,
933:                                                  PetscInt n,const PetscInt idx[],PetscInt *nout,PetscInt idxout[])
934: {

939:   if (!mapping->data) {
940:     ISGlobalToLocalMappingSetUp(mapping);
941:   }
942:   (*mapping->ops->globaltolocalmappingapplyblock)(mapping,type,n,idx,nout,idxout);
943:   return(0);
944: }


947: /*@C
948:     ISLocalToGlobalMappingGetBlockInfo - Gets the neighbor information for each processor and
949:      each index shared by more than one processor

951:     Collective on ISLocalToGlobalMapping

953:     Input Parameters:
954: .   mapping - the mapping from local to global indexing

956:     Output Parameter:
957: +   nproc - number of processors that are connected to this one
958: .   proc - neighboring processors
959: .   numproc - number of indices for each subdomain (processor)
960: -   indices - indices of nodes (in local numbering) shared with neighbors (sorted by global numbering)

962:     Level: advanced

964:     Concepts: mapping^local to global

966:     Fortran Usage:
967: $        ISLocalToGlobalMpngGetInfoSize(ISLocalToGlobalMapping,PetscInt nproc,PetscInt numprocmax,ierr) followed by
968: $        ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping,PetscInt nproc, PetscInt procs[nproc],PetscInt numprocs[nproc],
969:           PetscInt indices[nproc][numprocmax],ierr)
970:         There is no ISLocalToGlobalMappingRestoreInfo() in Fortran. You must make sure that procs[], numprocs[] and
971:         indices[][] are large enough arrays, either by allocating them dynamically or defining static ones large enough.


974: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
975:           ISLocalToGlobalMappingRestoreInfo()
976: @*/
977: PetscErrorCode  ISLocalToGlobalMappingGetBlockInfo(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
978: {

983:   if (mapping->info_cached) {
984:     *nproc    = mapping->info_nproc;
985:     *procs    = mapping->info_procs;
986:     *numprocs = mapping->info_numprocs;
987:     *indices  = mapping->info_indices;
988:   } else {
989:     ISLocalToGlobalMappingGetBlockInfo_Private(mapping,nproc,procs,numprocs,indices);
990:   }
991:   return(0);
992: }

994: static PetscErrorCode  ISLocalToGlobalMappingGetBlockInfo_Private(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
995: {
997:   PetscMPIInt    size,rank,tag1,tag2,tag3,*len,*source,imdex;
998:   PetscInt       i,n = mapping->n,Ng,ng,max = 0,*lindices = mapping->indices;
999:   PetscInt       *nprocs,*owner,nsends,*sends,j,*starts,nmax,nrecvs,*recvs,proc;
1000:   PetscInt       cnt,scale,*ownedsenders,*nownedsenders,rstart,nowned;
1001:   PetscInt       node,nownedm,nt,*sends2,nsends2,*starts2,*lens2,*dest,nrecvs2,*starts3,*recvs2,k,*bprocs,*tmp;
1002:   PetscInt       first_procs,first_numprocs,*first_indices;
1003:   MPI_Request    *recv_waits,*send_waits;
1004:   MPI_Status     recv_status,*send_status,*recv_statuses;
1005:   MPI_Comm       comm;
1006:   PetscBool      debug = PETSC_FALSE;

1009:   PetscObjectGetComm((PetscObject)mapping,&comm);
1010:   MPI_Comm_size(comm,&size);
1011:   MPI_Comm_rank(comm,&rank);
1012:   if (size == 1) {
1013:     *nproc         = 0;
1014:     *procs         = NULL;
1015:     PetscNew(numprocs);
1016:     (*numprocs)[0] = 0;
1017:     PetscNew(indices);
1018:     (*indices)[0]  = NULL;
1019:     /* save info for reuse */
1020:     mapping->info_nproc = *nproc;
1021:     mapping->info_procs = *procs;
1022:     mapping->info_numprocs = *numprocs;
1023:     mapping->info_indices = *indices;
1024:     mapping->info_cached = PETSC_TRUE;
1025:     return(0);
1026:   }

1028:   PetscOptionsGetBool(((PetscObject)mapping)->options,NULL,"-islocaltoglobalmappinggetinfo_debug",&debug,NULL);

1030:   /*
1031:     Notes on ISLocalToGlobalMappingGetBlockInfo

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

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

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

1041:     nontrivial locally owned node - node that is not in the interior (i.e. has more than one
1042:            local subdomain
1043:   */
1044:   PetscObjectGetNewTag((PetscObject)mapping,&tag1);
1045:   PetscObjectGetNewTag((PetscObject)mapping,&tag2);
1046:   PetscObjectGetNewTag((PetscObject)mapping,&tag3);

1048:   for (i=0; i<n; i++) {
1049:     if (lindices[i] > max) max = lindices[i];
1050:   }
1051:   MPIU_Allreduce(&max,&Ng,1,MPIU_INT,MPI_MAX,comm);
1052:   Ng++;
1053:   MPI_Comm_size(comm,&size);
1054:   MPI_Comm_rank(comm,&rank);
1055:   scale  = Ng/size + 1;
1056:   ng     = scale; if (rank == size-1) ng = Ng - scale*(size-1); ng = PetscMax(1,ng);
1057:   rstart = scale*rank;

1059:   /* determine ownership ranges of global indices */
1060:   PetscMalloc1(2*size,&nprocs);
1061:   PetscMemzero(nprocs,2*size*sizeof(PetscInt));

1063:   /* determine owners of each local node  */
1064:   PetscMalloc1(n,&owner);
1065:   for (i=0; i<n; i++) {
1066:     proc             = lindices[i]/scale; /* processor that globally owns this index */
1067:     nprocs[2*proc+1] = 1;                 /* processor globally owns at least one of ours */
1068:     owner[i]         = proc;
1069:     nprocs[2*proc]++;                     /* count of how many that processor globally owns of ours */
1070:   }
1071:   nsends = 0; for (i=0; i<size; i++) nsends += nprocs[2*i+1];
1072:   PetscInfo1(mapping,"Number of global owners for my local data %D\n",nsends);

1074:   /* inform other processors of number of messages and max length*/
1075:   PetscMaxSum(comm,nprocs,&nmax,&nrecvs);
1076:   PetscInfo1(mapping,"Number of local owners for my global data %D\n",nrecvs);

1078:   /* post receives for owned rows */
1079:   PetscMalloc1((2*nrecvs+1)*(nmax+1),&recvs);
1080:   PetscMalloc1(nrecvs+1,&recv_waits);
1081:   for (i=0; i<nrecvs; i++) {
1082:     MPI_Irecv(recvs+2*nmax*i,2*nmax,MPIU_INT,MPI_ANY_SOURCE,tag1,comm,recv_waits+i);
1083:   }

1085:   /* pack messages containing lists of local nodes to owners */
1086:   PetscMalloc1(2*n+1,&sends);
1087:   PetscMalloc1(size+1,&starts);
1088:   starts[0] = 0;
1089:   for (i=1; i<size; i++) starts[i] = starts[i-1] + 2*nprocs[2*i-2];
1090:   for (i=0; i<n; i++) {
1091:     sends[starts[owner[i]]++] = lindices[i];
1092:     sends[starts[owner[i]]++] = i;
1093:   }
1094:   PetscFree(owner);
1095:   starts[0] = 0;
1096:   for (i=1; i<size; i++) starts[i] = starts[i-1] + 2*nprocs[2*i-2];

1098:   /* send the messages */
1099:   PetscMalloc1(nsends+1,&send_waits);
1100:   PetscMalloc1(nsends+1,&dest);
1101:   cnt = 0;
1102:   for (i=0; i<size; i++) {
1103:     if (nprocs[2*i]) {
1104:       MPI_Isend(sends+starts[i],2*nprocs[2*i],MPIU_INT,i,tag1,comm,send_waits+cnt);
1105:       dest[cnt] = i;
1106:       cnt++;
1107:     }
1108:   }
1109:   PetscFree(starts);

1111:   /* wait on receives */
1112:   PetscMalloc1(nrecvs+1,&source);
1113:   PetscMalloc1(nrecvs+1,&len);
1114:   cnt  = nrecvs;
1115:   PetscMalloc1(ng+1,&nownedsenders);
1116:   PetscMemzero(nownedsenders,ng*sizeof(PetscInt));
1117:   while (cnt) {
1118:     MPI_Waitany(nrecvs,recv_waits,&imdex,&recv_status);
1119:     /* unpack receives into our local space */
1120:     MPI_Get_count(&recv_status,MPIU_INT,&len[imdex]);
1121:     source[imdex] = recv_status.MPI_SOURCE;
1122:     len[imdex]    = len[imdex]/2;
1123:     /* count how many local owners for each of my global owned indices */
1124:     for (i=0; i<len[imdex]; i++) nownedsenders[recvs[2*imdex*nmax+2*i]-rstart]++;
1125:     cnt--;
1126:   }
1127:   PetscFree(recv_waits);

1129:   /* count how many globally owned indices are on an edge multiplied by how many processors own them. */
1130:   nowned  = 0;
1131:   nownedm = 0;
1132:   for (i=0; i<ng; i++) {
1133:     if (nownedsenders[i] > 1) {nownedm += nownedsenders[i]; nowned++;}
1134:   }

1136:   /* create single array to contain rank of all local owners of each globally owned index */
1137:   PetscMalloc1(nownedm+1,&ownedsenders);
1138:   PetscMalloc1(ng+1,&starts);
1139:   starts[0] = 0;
1140:   for (i=1; i<ng; i++) {
1141:     if (nownedsenders[i-1] > 1) starts[i] = starts[i-1] + nownedsenders[i-1];
1142:     else starts[i] = starts[i-1];
1143:   }

1145:   /* for each nontrival globally owned node list all arriving processors */
1146:   for (i=0; i<nrecvs; i++) {
1147:     for (j=0; j<len[i]; j++) {
1148:       node = recvs[2*i*nmax+2*j]-rstart;
1149:       if (nownedsenders[node] > 1) ownedsenders[starts[node]++] = source[i];
1150:     }
1151:   }

1153:   if (debug) { /* -----------------------------------  */
1154:     starts[0] = 0;
1155:     for (i=1; i<ng; i++) {
1156:       if (nownedsenders[i-1] > 1) starts[i] = starts[i-1] + nownedsenders[i-1];
1157:       else starts[i] = starts[i-1];
1158:     }
1159:     for (i=0; i<ng; i++) {
1160:       if (nownedsenders[i] > 1) {
1161:         PetscSynchronizedPrintf(comm,"[%d] global node %D local owner processors: ",rank,i+rstart);
1162:         for (j=0; j<nownedsenders[i]; j++) {
1163:           PetscSynchronizedPrintf(comm,"%D ",ownedsenders[starts[i]+j]);
1164:         }
1165:         PetscSynchronizedPrintf(comm,"\n");
1166:       }
1167:     }
1168:     PetscSynchronizedFlush(comm,PETSC_STDOUT);
1169:   } /* -----------------------------------  */

1171:   /* wait on original sends */
1172:   if (nsends) {
1173:     PetscMalloc1(nsends,&send_status);
1174:     MPI_Waitall(nsends,send_waits,send_status);
1175:     PetscFree(send_status);
1176:   }
1177:   PetscFree(send_waits);
1178:   PetscFree(sends);
1179:   PetscFree(nprocs);

1181:   /* pack messages to send back to local owners */
1182:   starts[0] = 0;
1183:   for (i=1; i<ng; i++) {
1184:     if (nownedsenders[i-1] > 1) starts[i] = starts[i-1] + nownedsenders[i-1];
1185:     else starts[i] = starts[i-1];
1186:   }
1187:   nsends2 = nrecvs;
1188:   PetscMalloc1(nsends2+1,&nprocs); /* length of each message */
1189:   for (i=0; i<nrecvs; i++) {
1190:     nprocs[i] = 1;
1191:     for (j=0; j<len[i]; j++) {
1192:       node = recvs[2*i*nmax+2*j]-rstart;
1193:       if (nownedsenders[node] > 1) nprocs[i] += 2 + nownedsenders[node];
1194:     }
1195:   }
1196:   nt = 0;
1197:   for (i=0; i<nsends2; i++) nt += nprocs[i];

1199:   PetscMalloc1(nt+1,&sends2);
1200:   PetscMalloc1(nsends2+1,&starts2);

1202:   starts2[0] = 0;
1203:   for (i=1; i<nsends2; i++) starts2[i] = starts2[i-1] + nprocs[i-1];
1204:   /*
1205:      Each message is 1 + nprocs[i] long, and consists of
1206:        (0) the number of nodes being sent back
1207:        (1) the local node number,
1208:        (2) the number of processors sharing it,
1209:        (3) the processors sharing it
1210:   */
1211:   for (i=0; i<nsends2; i++) {
1212:     cnt = 1;
1213:     sends2[starts2[i]] = 0;
1214:     for (j=0; j<len[i]; j++) {
1215:       node = recvs[2*i*nmax+2*j]-rstart;
1216:       if (nownedsenders[node] > 1) {
1217:         sends2[starts2[i]]++;
1218:         sends2[starts2[i]+cnt++] = recvs[2*i*nmax+2*j+1];
1219:         sends2[starts2[i]+cnt++] = nownedsenders[node];
1220:         PetscMemcpy(&sends2[starts2[i]+cnt],&ownedsenders[starts[node]],nownedsenders[node]*sizeof(PetscInt));
1221:         cnt += nownedsenders[node];
1222:       }
1223:     }
1224:   }

1226:   /* receive the message lengths */
1227:   nrecvs2 = nsends;
1228:   PetscMalloc1(nrecvs2+1,&lens2);
1229:   PetscMalloc1(nrecvs2+1,&starts3);
1230:   PetscMalloc1(nrecvs2+1,&recv_waits);
1231:   for (i=0; i<nrecvs2; i++) {
1232:     MPI_Irecv(&lens2[i],1,MPIU_INT,dest[i],tag2,comm,recv_waits+i);
1233:   }

1235:   /* send the message lengths */
1236:   for (i=0; i<nsends2; i++) {
1237:     MPI_Send(&nprocs[i],1,MPIU_INT,source[i],tag2,comm);
1238:   }

1240:   /* wait on receives of lens */
1241:   if (nrecvs2) {
1242:     PetscMalloc1(nrecvs2,&recv_statuses);
1243:     MPI_Waitall(nrecvs2,recv_waits,recv_statuses);
1244:     PetscFree(recv_statuses);
1245:   }
1246:   PetscFree(recv_waits);

1248:   starts3[0] = 0;
1249:   nt         = 0;
1250:   for (i=0; i<nrecvs2-1; i++) {
1251:     starts3[i+1] = starts3[i] + lens2[i];
1252:     nt          += lens2[i];
1253:   }
1254:   if (nrecvs2) nt += lens2[nrecvs2-1];

1256:   PetscMalloc1(nt+1,&recvs2);
1257:   PetscMalloc1(nrecvs2+1,&recv_waits);
1258:   for (i=0; i<nrecvs2; i++) {
1259:     MPI_Irecv(recvs2+starts3[i],lens2[i],MPIU_INT,dest[i],tag3,comm,recv_waits+i);
1260:   }

1262:   /* send the messages */
1263:   PetscMalloc1(nsends2+1,&send_waits);
1264:   for (i=0; i<nsends2; i++) {
1265:     MPI_Isend(sends2+starts2[i],nprocs[i],MPIU_INT,source[i],tag3,comm,send_waits+i);
1266:   }

1268:   /* wait on receives */
1269:   if (nrecvs2) {
1270:     PetscMalloc1(nrecvs2,&recv_statuses);
1271:     MPI_Waitall(nrecvs2,recv_waits,recv_statuses);
1272:     PetscFree(recv_statuses);
1273:   }
1274:   PetscFree(recv_waits);
1275:   PetscFree(nprocs);

1277:   if (debug) { /* -----------------------------------  */
1278:     cnt = 0;
1279:     for (i=0; i<nrecvs2; i++) {
1280:       nt = recvs2[cnt++];
1281:       for (j=0; j<nt; j++) {
1282:         PetscSynchronizedPrintf(comm,"[%d] local node %D number of subdomains %D: ",rank,recvs2[cnt],recvs2[cnt+1]);
1283:         for (k=0; k<recvs2[cnt+1]; k++) {
1284:           PetscSynchronizedPrintf(comm,"%D ",recvs2[cnt+2+k]);
1285:         }
1286:         cnt += 2 + recvs2[cnt+1];
1287:         PetscSynchronizedPrintf(comm,"\n");
1288:       }
1289:     }
1290:     PetscSynchronizedFlush(comm,PETSC_STDOUT);
1291:   } /* -----------------------------------  */

1293:   /* count number subdomains for each local node */
1294:   PetscMalloc1(size,&nprocs);
1295:   PetscMemzero(nprocs,size*sizeof(PetscInt));
1296:   cnt  = 0;
1297:   for (i=0; i<nrecvs2; i++) {
1298:     nt = recvs2[cnt++];
1299:     for (j=0; j<nt; j++) {
1300:       for (k=0; k<recvs2[cnt+1]; k++) nprocs[recvs2[cnt+2+k]]++;
1301:       cnt += 2 + recvs2[cnt+1];
1302:     }
1303:   }
1304:   nt = 0; for (i=0; i<size; i++) nt += (nprocs[i] > 0);
1305:   *nproc    = nt;
1306:   PetscMalloc1(nt+1,procs);
1307:   PetscMalloc1(nt+1,numprocs);
1308:   PetscMalloc1(nt+1,indices);
1309:   for (i=0;i<nt+1;i++) (*indices)[i]=NULL;
1310:   PetscMalloc1(size,&bprocs);
1311:   cnt  = 0;
1312:   for (i=0; i<size; i++) {
1313:     if (nprocs[i] > 0) {
1314:       bprocs[i]        = cnt;
1315:       (*procs)[cnt]    = i;
1316:       (*numprocs)[cnt] = nprocs[i];
1317:       PetscMalloc1(nprocs[i],&(*indices)[cnt]);
1318:       cnt++;
1319:     }
1320:   }

1322:   /* make the list of subdomains for each nontrivial local node */
1323:   PetscMemzero(*numprocs,nt*sizeof(PetscInt));
1324:   cnt  = 0;
1325:   for (i=0; i<nrecvs2; i++) {
1326:     nt = recvs2[cnt++];
1327:     for (j=0; j<nt; j++) {
1328:       for (k=0; k<recvs2[cnt+1]; k++) (*indices)[bprocs[recvs2[cnt+2+k]]][(*numprocs)[bprocs[recvs2[cnt+2+k]]]++] = recvs2[cnt];
1329:       cnt += 2 + recvs2[cnt+1];
1330:     }
1331:   }
1332:   PetscFree(bprocs);
1333:   PetscFree(recvs2);

1335:   /* sort the node indexing by their global numbers */
1336:   nt = *nproc;
1337:   for (i=0; i<nt; i++) {
1338:     PetscMalloc1((*numprocs)[i],&tmp);
1339:     for (j=0; j<(*numprocs)[i]; j++) tmp[j] = lindices[(*indices)[i][j]];
1340:     PetscSortIntWithArray((*numprocs)[i],tmp,(*indices)[i]);
1341:     PetscFree(tmp);
1342:   }

1344:   if (debug) { /* -----------------------------------  */
1345:     nt = *nproc;
1346:     for (i=0; i<nt; i++) {
1347:       PetscSynchronizedPrintf(comm,"[%d] subdomain %D number of indices %D: ",rank,(*procs)[i],(*numprocs)[i]);
1348:       for (j=0; j<(*numprocs)[i]; j++) {
1349:         PetscSynchronizedPrintf(comm,"%D ",(*indices)[i][j]);
1350:       }
1351:       PetscSynchronizedPrintf(comm,"\n");
1352:     }
1353:     PetscSynchronizedFlush(comm,PETSC_STDOUT);
1354:   } /* -----------------------------------  */

1356:   /* wait on sends */
1357:   if (nsends2) {
1358:     PetscMalloc1(nsends2,&send_status);
1359:     MPI_Waitall(nsends2,send_waits,send_status);
1360:     PetscFree(send_status);
1361:   }

1363:   PetscFree(starts3);
1364:   PetscFree(dest);
1365:   PetscFree(send_waits);

1367:   PetscFree(nownedsenders);
1368:   PetscFree(ownedsenders);
1369:   PetscFree(starts);
1370:   PetscFree(starts2);
1371:   PetscFree(lens2);

1373:   PetscFree(source);
1374:   PetscFree(len);
1375:   PetscFree(recvs);
1376:   PetscFree(nprocs);
1377:   PetscFree(sends2);

1379:   /* put the information about myself as the first entry in the list */
1380:   first_procs    = (*procs)[0];
1381:   first_numprocs = (*numprocs)[0];
1382:   first_indices  = (*indices)[0];
1383:   for (i=0; i<*nproc; i++) {
1384:     if ((*procs)[i] == rank) {
1385:       (*procs)[0]    = (*procs)[i];
1386:       (*numprocs)[0] = (*numprocs)[i];
1387:       (*indices)[0]  = (*indices)[i];
1388:       (*procs)[i]    = first_procs;
1389:       (*numprocs)[i] = first_numprocs;
1390:       (*indices)[i]  = first_indices;
1391:       break;
1392:     }
1393:   }

1395:   /* save info for reuse */
1396:   mapping->info_nproc = *nproc;
1397:   mapping->info_procs = *procs;
1398:   mapping->info_numprocs = *numprocs;
1399:   mapping->info_indices = *indices;
1400:   mapping->info_cached = PETSC_TRUE;
1401:   return(0);
1402: }

1404: /*@C
1405:     ISLocalToGlobalMappingRestoreBlockInfo - Frees the memory allocated by ISLocalToGlobalMappingGetBlockInfo()

1407:     Collective on ISLocalToGlobalMapping

1409:     Input Parameters:
1410: .   mapping - the mapping from local to global indexing

1412:     Output Parameter:
1413: +   nproc - number of processors that are connected to this one
1414: .   proc - neighboring processors
1415: .   numproc - number of indices for each processor
1416: -   indices - indices of local nodes shared with neighbor (sorted by global numbering)

1418:     Level: advanced

1420: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
1421:           ISLocalToGlobalMappingGetInfo()
1422: @*/
1423: PetscErrorCode  ISLocalToGlobalMappingRestoreBlockInfo(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
1424: {

1429:   if (mapping->info_free) {
1430:     PetscFree(*numprocs);
1431:     if (*indices) {
1432:       PetscInt i;

1434:       PetscFree((*indices)[0]);
1435:       for (i=1; i<*nproc; i++) {
1436:         PetscFree((*indices)[i]);
1437:       }
1438:       PetscFree(*indices);
1439:     }
1440:   }
1441:   *nproc    = 0;
1442:   *procs    = NULL;
1443:   *numprocs = NULL;
1444:   *indices  = NULL;
1445:   return(0);
1446: }

1448: /*@C
1449:     ISLocalToGlobalMappingGetInfo - Gets the neighbor information for each processor and
1450:      each index shared by more than one processor

1452:     Collective on ISLocalToGlobalMapping

1454:     Input Parameters:
1455: .   mapping - the mapping from local to global indexing

1457:     Output Parameter:
1458: +   nproc - number of processors that are connected to this one
1459: .   proc - neighboring processors
1460: .   numproc - number of indices for each subdomain (processor)
1461: -   indices - indices of nodes (in local numbering) shared with neighbors (sorted by global numbering)

1463:     Level: advanced

1465:     Concepts: mapping^local to global

1467:     Notes: The user needs to call ISLocalToGlobalMappingRestoreInfo when the data is no longer needed.

1469:     Fortran Usage:
1470: $        ISLocalToGlobalMpngGetInfoSize(ISLocalToGlobalMapping,PetscInt nproc,PetscInt numprocmax,ierr) followed by
1471: $        ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping,PetscInt nproc, PetscInt procs[nproc],PetscInt numprocs[nproc],
1472:           PetscInt indices[nproc][numprocmax],ierr)
1473:         There is no ISLocalToGlobalMappingRestoreInfo() in Fortran. You must make sure that procs[], numprocs[] and
1474:         indices[][] are large enough arrays, either by allocating them dynamically or defining static ones large enough.


1477: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
1478:           ISLocalToGlobalMappingRestoreInfo()
1479: @*/
1480: PetscErrorCode  ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
1481: {
1483:   PetscInt       **bindices = NULL,*bnumprocs = NULL,bs = mapping->bs,i,j,k;

1487:   ISLocalToGlobalMappingGetBlockInfo(mapping,nproc,procs,&bnumprocs,&bindices);
1488:   if (bs > 1) { /* we need to expand the cached info */
1489:     PetscCalloc1(*nproc,&*indices);
1490:     PetscCalloc1(*nproc,&*numprocs);
1491:     for (i=0; i<*nproc; i++) {
1492:       PetscMalloc1(bs*bnumprocs[i],&(*indices)[i]);
1493:       for (j=0; j<bnumprocs[i]; j++) {
1494:         for (k=0; k<bs; k++) {
1495:           (*indices)[i][j*bs+k] = bs*bindices[i][j] + k;
1496:         }
1497:       }
1498:       (*numprocs)[i] = bnumprocs[i]*bs;
1499:     }
1500:     mapping->info_free = PETSC_TRUE;
1501:   } else {
1502:     *numprocs = bnumprocs;
1503:     *indices  = bindices;
1504:   }
1505:   return(0);
1506: }

1508: /*@C
1509:     ISLocalToGlobalMappingRestoreInfo - Frees the memory allocated by ISLocalToGlobalMappingGetInfo()

1511:     Collective on ISLocalToGlobalMapping

1513:     Input Parameters:
1514: .   mapping - the mapping from local to global indexing

1516:     Output Parameter:
1517: +   nproc - number of processors that are connected to this one
1518: .   proc - neighboring processors
1519: .   numproc - number of indices for each processor
1520: -   indices - indices of local nodes shared with neighbor (sorted by global numbering)

1522:     Level: advanced

1524: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
1525:           ISLocalToGlobalMappingGetInfo()
1526: @*/
1527: PetscErrorCode  ISLocalToGlobalMappingRestoreInfo(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
1528: {

1532:   ISLocalToGlobalMappingRestoreBlockInfo(mapping,nproc,procs,numprocs,indices);
1533:   return(0);
1534: }

1536: /*@C
1537:     ISLocalToGlobalMappingGetNodeInfo - Gets the neighbor information for each node

1539:     Collective on ISLocalToGlobalMapping

1541:     Input Parameters:
1542: .   mapping - the mapping from local to global indexing

1544:     Output Parameter:
1545: +   nnodes - number of local nodes (same ISLocalToGlobalMappingGetSize())
1546: .   count - number of neighboring processors per node
1547: -   indices - indices of processes sharing the node (sorted)

1549:     Level: advanced

1551:     Concepts: mapping^local to global

1553:     Notes: The user needs to call ISLocalToGlobalMappingRestoreInfo when the data is no longer needed.

1555: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
1556:           ISLocalToGlobalMappingGetInfo(), ISLocalToGlobalMappingRestoreNodeInfo()
1557: @*/
1558: PetscErrorCode  ISLocalToGlobalMappingGetNodeInfo(ISLocalToGlobalMapping mapping,PetscInt *nnodes,PetscInt *count[],PetscInt **indices[])
1559: {
1560:   PetscInt       n;

1565:   ISLocalToGlobalMappingGetSize(mapping,&n);
1566:   if (!mapping->info_nodec) {
1567:     PetscInt i,m,n_neigh,*neigh,*n_shared,**shared;

1569:     PetscCalloc1(n+1,&mapping->info_nodec); /* always allocate to flag setup */
1570:     PetscMalloc1(n,&mapping->info_nodei);
1571:     ISLocalToGlobalMappingGetInfo(mapping,&n_neigh,&neigh,&n_shared,&shared);
1572:     for (i=0,m=0;i<n;i++) { mapping->info_nodec[i] = 1; m++; }
1573:     for (i=1;i<n_neigh;i++) {
1574:       PetscInt j;

1576:       m += n_shared[i];
1577:       for (j=0;j<n_shared[i];j++) mapping->info_nodec[shared[i][j]] += 1;
1578:     }
1579:     if (n) { PetscMalloc1(m,&mapping->info_nodei[0]); }
1580:     for (i=1;i<n;i++) mapping->info_nodei[i] = mapping->info_nodei[i-1] + mapping->info_nodec[i-1];
1581:     PetscMemzero(mapping->info_nodec,n*sizeof(PetscInt));
1582:     for (i=0;i<n;i++) { mapping->info_nodec[i] = 1; mapping->info_nodei[i][0] = neigh[0]; }
1583:     for (i=1;i<n_neigh;i++) {
1584:       PetscInt j;

1586:       for (j=0;j<n_shared[i];j++) {
1587:         PetscInt k = shared[i][j];

1589:         mapping->info_nodei[k][mapping->info_nodec[k]] = neigh[i];
1590:         mapping->info_nodec[k] += 1;
1591:       }
1592:     }
1593:     for (i=0;i<n;i++) { PetscSortRemoveDupsInt(&mapping->info_nodec[i],mapping->info_nodei[i]); }
1594:     ISLocalToGlobalMappingRestoreInfo(mapping,&n_neigh,&neigh,&n_shared,&shared);
1595:   }
1596:   if (nnodes)  *nnodes  = n;
1597:   if (count)   *count   = mapping->info_nodec;
1598:   if (indices) *indices = mapping->info_nodei;
1599:   return(0);
1600: }

1602: /*@C
1603:     ISLocalToGlobalMappingRestoreNodeInfo - Frees the memory allocated by ISLocalToGlobalMappingGetNodeInfo()

1605:     Collective on ISLocalToGlobalMapping

1607:     Input Parameters:
1608: .   mapping - the mapping from local to global indexing

1610:     Output Parameter:
1611: +   nnodes - number of local nodes
1612: .   count - number of neighboring processors per node
1613: -   indices - indices of processes sharing the node (sorted)

1615:     Level: advanced

1617: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
1618:           ISLocalToGlobalMappingGetInfo()
1619: @*/
1620: PetscErrorCode  ISLocalToGlobalMappingRestoreNodeInfo(ISLocalToGlobalMapping mapping,PetscInt *nnodes,PetscInt *count[],PetscInt **indices[])
1621: {
1624:   if (nnodes)  *nnodes  = 0;
1625:   if (count)   *count   = NULL;
1626:   if (indices) *indices = NULL;
1627:   return(0);
1628: }

1630: /*@C
1631:    ISLocalToGlobalMappingGetIndices - Get global indices for every local point that is mapped

1633:    Not Collective

1635:    Input Arguments:
1636: . ltog - local to global mapping

1638:    Output Arguments:
1639: . array - array of indices, the length of this array may be obtained with ISLocalToGlobalMappingGetSize()

1641:    Level: advanced

1643:    Notes:
1644:     ISLocalToGlobalMappingGetSize() returns the length the this array

1646: .seealso: ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingRestoreIndices(), ISLocalToGlobalMappingGetBlockIndices(), ISLocalToGlobalMappingRestoreBlockIndices()
1647: @*/
1648: PetscErrorCode  ISLocalToGlobalMappingGetIndices(ISLocalToGlobalMapping ltog,const PetscInt **array)
1649: {
1653:   if (ltog->bs == 1) {
1654:     *array = ltog->indices;
1655:   } else {
1656:     PetscInt       *jj,k,i,j,n = ltog->n, bs = ltog->bs;
1657:     const PetscInt *ii;

1660:     PetscMalloc1(bs*n,&jj);
1661:     *array = jj;
1662:     k    = 0;
1663:     ii   = ltog->indices;
1664:     for (i=0; i<n; i++)
1665:       for (j=0; j<bs; j++)
1666:         jj[k++] = bs*ii[i] + j;
1667:   }
1668:   return(0);
1669: }

1671: /*@C
1672:    ISLocalToGlobalMappingRestoreIndices - Restore indices obtained with ISLocalToGlobalMappingGetIndices()

1674:    Not Collective

1676:    Input Arguments:
1677: + ltog - local to global mapping
1678: - array - array of indices

1680:    Level: advanced

1682: .seealso: ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingGetIndices()
1683: @*/
1684: PetscErrorCode  ISLocalToGlobalMappingRestoreIndices(ISLocalToGlobalMapping ltog,const PetscInt **array)
1685: {
1689:   if (ltog->bs == 1 && *array != ltog->indices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_BADPTR,"Trying to return mismatched pointer");

1691:   if (ltog->bs > 1) {
1693:     PetscFree(*(void**)array);
1694:   }
1695:   return(0);
1696: }

1698: /*@C
1699:    ISLocalToGlobalMappingGetBlockIndices - Get global indices for every local block

1701:    Not Collective

1703:    Input Arguments:
1704: . ltog - local to global mapping

1706:    Output Arguments:
1707: . array - array of indices

1709:    Level: advanced

1711: .seealso: ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingRestoreBlockIndices()
1712: @*/
1713: PetscErrorCode  ISLocalToGlobalMappingGetBlockIndices(ISLocalToGlobalMapping ltog,const PetscInt **array)
1714: {
1718:   *array = ltog->indices;
1719:   return(0);
1720: }

1722: /*@C
1723:    ISLocalToGlobalMappingRestoreBlockIndices - Restore indices obtained with ISLocalToGlobalMappingGetBlockIndices()

1725:    Not Collective

1727:    Input Arguments:
1728: + ltog - local to global mapping
1729: - array - array of indices

1731:    Level: advanced

1733: .seealso: ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingGetIndices()
1734: @*/
1735: PetscErrorCode  ISLocalToGlobalMappingRestoreBlockIndices(ISLocalToGlobalMapping ltog,const PetscInt **array)
1736: {
1740:   if (*array != ltog->indices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_BADPTR,"Trying to return mismatched pointer");
1741:   *array = NULL;
1742:   return(0);
1743: }

1745: /*@C
1746:    ISLocalToGlobalMappingConcatenate - Create a new mapping that concatenates a list of mappings

1748:    Not Collective

1750:    Input Arguments:
1751: + comm - communicator for the new mapping, must contain the communicator of every mapping to concatenate
1752: . n - number of mappings to concatenate
1753: - ltogs - local to global mappings

1755:    Output Arguments:
1756: . ltogcat - new mapping

1758:    Note: this currently always returns a mapping with block size of 1

1760:    Developer Note: If all the input mapping have the same block size we could easily handle that as a special case

1762:    Level: advanced

1764: .seealso: ISLocalToGlobalMappingCreate()
1765: @*/
1766: PetscErrorCode ISLocalToGlobalMappingConcatenate(MPI_Comm comm,PetscInt n,const ISLocalToGlobalMapping ltogs[],ISLocalToGlobalMapping *ltogcat)
1767: {
1768:   PetscInt       i,cnt,m,*idx;

1772:   if (n < 0) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have a non-negative number of mappings, given %D",n);
1776:   for (cnt=0,i=0; i<n; i++) {
1777:     ISLocalToGlobalMappingGetSize(ltogs[i],&m);
1778:     cnt += m;
1779:   }
1780:   PetscMalloc1(cnt,&idx);
1781:   for (cnt=0,i=0; i<n; i++) {
1782:     const PetscInt *subidx;
1783:     ISLocalToGlobalMappingGetSize(ltogs[i],&m);
1784:     ISLocalToGlobalMappingGetIndices(ltogs[i],&subidx);
1785:     PetscMemcpy(&idx[cnt],subidx,m*sizeof(PetscInt));
1786:     ISLocalToGlobalMappingRestoreIndices(ltogs[i],&subidx);
1787:     cnt += m;
1788:   }
1789:   ISLocalToGlobalMappingCreate(comm,1,cnt,idx,PETSC_OWN_POINTER,ltogcat);
1790:   return(0);
1791: }

1793: /*MC
1794:       ISLOCALTOGLOBALMAPPINGBASIC - basic implementation of the ISLocalToGlobalMapping object. When ISGlobalToLocalMappingApply() is
1795:                                     used this is good for only small and moderate size problems.

1797:    Options Database Keys:
1798: +   -islocaltoglobalmapping_type basic - select this method

1800:    Level: beginner

1802: .seealso:  ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingSetType(), ISLOCALTOGLOBALMAPPINGHASH
1803: M*/
1804: PETSC_EXTERN PetscErrorCode ISLocalToGlobalMappingCreate_Basic(ISLocalToGlobalMapping ltog)
1805: {
1807:   ltog->ops->globaltolocalmappingapply      = ISGlobalToLocalMappingApply_Basic;
1808:   ltog->ops->globaltolocalmappingsetup      = ISGlobalToLocalMappingSetUp_Basic;
1809:   ltog->ops->globaltolocalmappingapplyblock = ISGlobalToLocalMappingApplyBlock_Basic;
1810:   ltog->ops->destroy                        = ISLocalToGlobalMappingDestroy_Basic;
1811:   return(0);
1812: }

1814: /*MC
1815:       ISLOCALTOGLOBALMAPPINGHASH - hash implementation of the ISLocalToGlobalMapping object. When ISGlobalToLocalMappingApply() is
1816:                                     used this is good for large memory problems.

1818:    Options Database Keys:
1819: +   -islocaltoglobalmapping_type hash - select this method


1822:    Notes:
1823:     This is selected automatically for large problems if the user does not set the type.

1825:    Level: beginner

1827: .seealso:  ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingSetType(), ISLOCALTOGLOBALMAPPINGHASH
1828: M*/
1829: PETSC_EXTERN PetscErrorCode ISLocalToGlobalMappingCreate_Hash(ISLocalToGlobalMapping ltog)
1830: {
1832:   ltog->ops->globaltolocalmappingapply      = ISGlobalToLocalMappingApply_Hash;
1833:   ltog->ops->globaltolocalmappingsetup      = ISGlobalToLocalMappingSetUp_Hash;
1834:   ltog->ops->globaltolocalmappingapplyblock = ISGlobalToLocalMappingApplyBlock_Hash;
1835:   ltog->ops->destroy                        = ISLocalToGlobalMappingDestroy_Hash;
1836:   return(0);
1837: }


1840: /*@C
1841:     ISLocalToGlobalMappingRegister -  Adds a method for applying a global to local mapping with an ISLocalToGlobalMapping

1843:    Not Collective

1845:    Input Parameters:
1846: +  sname - name of a new method
1847: -  routine_create - routine to create method context

1849:    Notes:
1850:    ISLocalToGlobalMappingRegister() may be called multiple times to add several user-defined mappings.

1852:    Sample usage:
1853: .vb
1854:    ISLocalToGlobalMappingRegister("my_mapper",MyCreate);
1855: .ve

1857:    Then, your mapping can be chosen with the procedural interface via
1858: $     ISLocalToGlobalMappingSetType(ltog,"my_mapper")
1859:    or at runtime via the option
1860: $     -islocaltoglobalmapping_type my_mapper

1862:    Level: advanced

1864: .keywords: ISLocalToGlobalMappingType, register

1866: .seealso: ISLocalToGlobalMappingRegisterAll(), ISLocalToGlobalMappingRegisterDestroy(), ISLOCALTOGLOBALMAPPINGBASIC, ISLOCALTOGLOBALMAPPINGHASH

1868: @*/
1869: PetscErrorCode  ISLocalToGlobalMappingRegister(const char sname[],PetscErrorCode (*function)(ISLocalToGlobalMapping))
1870: {

1874:   ISInitializePackage();
1875:   PetscFunctionListAdd(&ISLocalToGlobalMappingList,sname,function);
1876:   return(0);
1877: }

1879: /*@C
1880:    ISLocalToGlobalMappingSetType - Builds ISLocalToGlobalMapping for a particular global to local mapping approach.

1882:    Logically Collective on ISLocalToGlobalMapping

1884:    Input Parameters:
1885: +  ltog - the ISLocalToGlobalMapping object
1886: -  type - a known method

1888:    Options Database Key:
1889: .  -islocaltoglobalmapping_type  <method> - Sets the method; use -help for a list
1890:     of available methods (for instance, basic or hash)

1892:    Notes:
1893:    See "petsc/include/petscis.h" for available methods

1895:   Normally, it is best to use the ISLocalToGlobalMappingSetFromOptions() command and
1896:   then set the ISLocalToGlobalMapping type from the options database rather than by using
1897:   this routine.

1899:   Level: intermediate

1901:   Developer Note: ISLocalToGlobalMappingRegister() is used to add new types to ISLocalToGlobalMappingList from which they
1902:   are accessed by ISLocalToGlobalMappingSetType().

1904: .keywords: ISLocalToGlobalMapping, set, method

1906: .seealso: PCSetType(), ISLocalToGlobalMappingType, ISLocalToGlobalMappingRegister(), ISLocalToGlobalMappingCreate()

1908: @*/
1909: PetscErrorCode  ISLocalToGlobalMappingSetType(ISLocalToGlobalMapping ltog, ISLocalToGlobalMappingType type)
1910: {
1911:   PetscErrorCode ierr,(*r)(ISLocalToGlobalMapping);
1912:   PetscBool      match;


1918:   PetscObjectTypeCompare((PetscObject)ltog,type,&match);
1919:   if (match) return(0);

1921:    PetscFunctionListFind(ISLocalToGlobalMappingList,type,&r);
1922:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)ltog),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested ISLocalToGlobalMapping type %s",type);
1923:   /* Destroy the previous private LTOG context */
1924:   if (ltog->ops->destroy) {
1925:     (*ltog->ops->destroy)(ltog);
1926:     ltog->ops->destroy = NULL;
1927:   }
1928:   PetscObjectChangeTypeName((PetscObject)ltog,type);
1929:   (*r)(ltog);
1930:   return(0);
1931: }

1933: PetscBool ISLocalToGlobalMappingRegisterAllCalled = PETSC_FALSE;

1935: /*@C
1936:   ISLocalToGlobalMappingRegisterAll - Registers all of the local to global mapping components in the IS package.

1938:   Not Collective

1940:   Level: advanced

1942: .keywords: IS, register, all
1943: .seealso:  ISRegister(),  ISLocalToGlobalRegister()
1944: @*/
1945: PetscErrorCode  ISLocalToGlobalMappingRegisterAll(void)
1946: {

1950:   if (ISLocalToGlobalMappingRegisterAllCalled) return(0);
1951:   ISLocalToGlobalMappingRegisterAllCalled = PETSC_TRUE;

1953:   ISLocalToGlobalMappingRegister(ISLOCALTOGLOBALMAPPINGBASIC, ISLocalToGlobalMappingCreate_Basic);
1954:   ISLocalToGlobalMappingRegister(ISLOCALTOGLOBALMAPPINGHASH, ISLocalToGlobalMappingCreate_Hash);
1955:   return(0);
1956: }