Actual source code: isltog.c

petsc-3.12.5 2020-03-29
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:     if (local >= 0) 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:     if (local >= 0) 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: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate()
213: @*/
214: PetscErrorCode  ISLocalToGlobalMappingDuplicate(ISLocalToGlobalMapping ltog,ISLocalToGlobalMapping* nltog)
215: {

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

224: /*@
225:     ISLocalToGlobalMappingGetSize - Gets the local size of a local to global mapping

227:     Not Collective

229:     Input Parameter:
230: .   ltog - local to global mapping

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

235:     Level: advanced

237: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate()
238: @*/
239: PetscErrorCode  ISLocalToGlobalMappingGetSize(ISLocalToGlobalMapping mapping,PetscInt *n)
240: {
244:   *n = mapping->bs*mapping->n;
245:   return(0);
246: }

248: /*@C
249:     ISLocalToGlobalMappingView - View a local to global mapping

251:     Not Collective

253:     Input Parameters:
254: +   ltog - local to global mapping
255: -   viewer - viewer

257:     Level: advanced

259: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate()
260: @*/
261: PetscErrorCode  ISLocalToGlobalMappingView(ISLocalToGlobalMapping mapping,PetscViewer viewer)
262: {
263:   PetscInt       i;
264:   PetscMPIInt    rank;
265:   PetscBool      iascii;

270:   if (!viewer) {
271:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mapping),&viewer);
272:   }

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

289: /*@
290:     ISLocalToGlobalMappingCreateIS - Creates a mapping between a local (0 to n)
291:     ordering and a global parallel ordering.

293:     Not collective

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

298:     Output Parameter:
299: .   mapping - new mapping data structure

301:     Notes:
302:     the block size of the IS determines the block size of the mapping
303:     Level: advanced

305: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingSetFromOptions()
306: @*/
307: PetscErrorCode  ISLocalToGlobalMappingCreateIS(IS is,ISLocalToGlobalMapping *mapping)
308: {
310:   PetscInt       n,bs;
311:   const PetscInt *indices;
312:   MPI_Comm       comm;
313:   PetscBool      isblock;


319:   PetscObjectGetComm((PetscObject)is,&comm);
320:   ISGetLocalSize(is,&n);
321:   PetscObjectTypeCompare((PetscObject)is,ISBLOCK,&isblock);
322:   if (!isblock) {
323:     ISGetIndices(is,&indices);
324:     ISLocalToGlobalMappingCreate(comm,1,n,indices,PETSC_COPY_VALUES,mapping);
325:     ISRestoreIndices(is,&indices);
326:   } else {
327:     ISGetBlockSize(is,&bs);
328:     ISBlockGetIndices(is,&indices);
329:     ISLocalToGlobalMappingCreate(comm,bs,n/bs,indices,PETSC_COPY_VALUES,mapping);
330:     ISBlockRestoreIndices(is,&indices);
331:   }
332:   return(0);
333: }

335: /*@C
336:     ISLocalToGlobalMappingCreateSF - Creates a mapping between a local (0 to n)
337:     ordering and a global parallel ordering.

339:     Collective

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

345:     Output Parameter:
346: .   mapping - new mapping data structure

348:     Level: advanced

350: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingSetFromOptions()
351: @*/
352: PetscErrorCode ISLocalToGlobalMappingCreateSF(PetscSF sf,PetscInt start,ISLocalToGlobalMapping *mapping)
353: {
355:   PetscInt       i,maxlocal,nroots,nleaves,*globals,*ltog;
356:   const PetscInt *ilocal;
357:   MPI_Comm       comm;


363:   PetscObjectGetComm((PetscObject)sf,&comm);
364:   PetscSFGetGraph(sf,&nroots,&nleaves,&ilocal,NULL);
365:   if (ilocal) {
366:     for (i=0,maxlocal=0; i<nleaves; i++) maxlocal = PetscMax(maxlocal,ilocal[i]+1);
367:   }
368:   else maxlocal = nleaves;
369:   PetscMalloc1(nroots,&globals);
370:   PetscMalloc1(maxlocal,&ltog);
371:   for (i=0; i<nroots; i++) globals[i] = start + i;
372:   for (i=0; i<maxlocal; i++) ltog[i] = -1;
373:   PetscSFBcastBegin(sf,MPIU_INT,globals,ltog);
374:   PetscSFBcastEnd(sf,MPIU_INT,globals,ltog);
375:   ISLocalToGlobalMappingCreate(comm,1,maxlocal,ltog,PETSC_OWN_POINTER,mapping);
376:   PetscFree(globals);
377:   return(0);
378: }

380: /*@
381:     ISLocalToGlobalMappingSetBlockSize - Sets the blocksize of the mapping

383:     Not collective

385:     Input Parameters:
386: +   mapping - mapping data structure
387: -   bs - the blocksize

389:     Level: advanced

391: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS()
392: @*/
393: PetscErrorCode  ISLocalToGlobalMappingSetBlockSize(ISLocalToGlobalMapping mapping,PetscInt bs)
394: {
395:   PetscInt       *nid;
396:   const PetscInt *oid;
397:   PetscInt       i,cn,on,obs,nn;

402:   if (bs < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid block size %D",bs);
403:   if (bs == mapping->bs) return(0);
404:   on  = mapping->n;
405:   obs = mapping->bs;
406:   oid = mapping->indices;
407:   nn  = (on*obs)/bs;
408:   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);

410:   PetscMalloc1(nn,&nid);
411:   ISLocalToGlobalMappingGetIndices(mapping,&oid);
412:   for (i=0;i<nn;i++) {
413:     PetscInt j;
414:     for (j=0,cn=0;j<bs-1;j++) {
415:       if (oid[i*bs+j] < 0) { cn++; continue; }
416:       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]);
417:     }
418:     if (oid[i*bs+j] < 0) cn++;
419:     if (cn) {
420:       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);
421:       nid[i] = -1;
422:     } else {
423:       nid[i] = oid[i*bs]/bs;
424:     }
425:   }
426:   ISLocalToGlobalMappingRestoreIndices(mapping,&oid);

428:   mapping->n           = nn;
429:   mapping->bs          = bs;
430:   PetscFree(mapping->indices);
431:   mapping->indices     = nid;
432:   mapping->globalstart = 0;
433:   mapping->globalend   = 0;

435:   /* reset the cached information */
436:   PetscFree(mapping->info_procs);
437:   PetscFree(mapping->info_numprocs);
438:   if (mapping->info_indices) {
439:     PetscInt i;

441:     PetscFree((mapping->info_indices)[0]);
442:     for (i=1; i<mapping->info_nproc; i++) {
443:       PetscFree(mapping->info_indices[i]);
444:     }
445:     PetscFree(mapping->info_indices);
446:   }
447:   mapping->info_cached = PETSC_FALSE;

449:   if (mapping->ops->destroy) {
450:     (*mapping->ops->destroy)(mapping);
451:   }
452:   return(0);
453: }

455: /*@
456:     ISLocalToGlobalMappingGetBlockSize - Gets the blocksize of the mapping
457:     ordering and a global parallel ordering.

459:     Not Collective

461:     Input Parameters:
462: .   mapping - mapping data structure

464:     Output Parameter:
465: .   bs - the blocksize

467:     Level: advanced

469: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS()
470: @*/
471: PetscErrorCode  ISLocalToGlobalMappingGetBlockSize(ISLocalToGlobalMapping mapping,PetscInt *bs)
472: {
475:   *bs = mapping->bs;
476:   return(0);
477: }

479: /*@
480:     ISLocalToGlobalMappingCreate - Creates a mapping between a local (0 to n)
481:     ordering and a global parallel ordering.

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

485:     Input Parameters:
486: +   comm - MPI communicator
487: .   bs - the block size
488: .   n - the number of local elements divided by the block size, or equivalently the number of block indices
489: .   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
490: -   mode - see PetscCopyMode

492:     Output Parameter:
493: .   mapping - new mapping data structure

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

498:     For "small" problems when using ISGlobalToLocalMappingApply() and ISGlobalToLocalMappingApplyBlock(), the ISLocalToGlobalMappingType of ISLOCALTOGLOBALMAPPINGBASIC will be used;
499:     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.
500:     Use ISLocalToGlobalMappingSetType() or call ISLocalToGlobalMappingSetFromOptions() with the option -islocaltoglobalmapping_type <basic,hash> to control which is used.

502:     Level: advanced

504: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingSetFromOptions(), ISLOCALTOGLOBALMAPPINGBASIC, ISLOCALTOGLOBALMAPPINGHASH
505:           ISLocalToGlobalMappingSetType(), ISLocalToGlobalMappingType
506: @*/
507: PetscErrorCode  ISLocalToGlobalMappingCreate(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscInt indices[],PetscCopyMode mode,ISLocalToGlobalMapping *mapping)
508: {
510:   PetscInt       *in;


516:   *mapping = NULL;
517:   ISInitializePackage();

519:   PetscHeaderCreate(*mapping,IS_LTOGM_CLASSID,"ISLocalToGlobalMapping","Local to global mapping","IS",
520:                            comm,ISLocalToGlobalMappingDestroy,ISLocalToGlobalMappingView);
521:   (*mapping)->n             = n;
522:   (*mapping)->bs            = bs;
523:   (*mapping)->info_cached   = PETSC_FALSE;
524:   (*mapping)->info_free     = PETSC_FALSE;
525:   (*mapping)->info_procs    = NULL;
526:   (*mapping)->info_numprocs = NULL;
527:   (*mapping)->info_indices  = NULL;
528:   (*mapping)->info_nodec    = NULL;
529:   (*mapping)->info_nodei    = NULL;

531:   (*mapping)->ops->globaltolocalmappingapply      = NULL;
532:   (*mapping)->ops->globaltolocalmappingapplyblock = NULL;
533:   (*mapping)->ops->destroy                        = NULL;
534:   if (mode == PETSC_COPY_VALUES) {
535:     PetscMalloc1(n,&in);
536:     PetscArraycpy(in,indices,n);
537:     (*mapping)->indices = in;
538:     PetscLogObjectMemory((PetscObject)*mapping,n*sizeof(PetscInt));
539:   } else if (mode == PETSC_OWN_POINTER) {
540:     (*mapping)->indices = (PetscInt*)indices;
541:     PetscLogObjectMemory((PetscObject)*mapping,n*sizeof(PetscInt));
542:   }
543:   else SETERRQ(comm,PETSC_ERR_SUP,"Cannot currently use PETSC_USE_POINTER");
544:   return(0);
545: }

547: PetscFunctionList ISLocalToGlobalMappingList = NULL;


550: /*@
551:    ISLocalToGlobalMappingSetFromOptions - Set mapping options from the options database.

553:    Not collective

555:    Input Parameters:
556: .  mapping - mapping data structure

558:    Level: advanced

560: @*/
561: PetscErrorCode ISLocalToGlobalMappingSetFromOptions(ISLocalToGlobalMapping mapping)
562: {
563:   PetscErrorCode             ierr;
564:   char                       type[256];
565:   ISLocalToGlobalMappingType defaulttype = "Not set";
566:   PetscBool                  flg;

570:   ISLocalToGlobalMappingRegisterAll();
571:   PetscObjectOptionsBegin((PetscObject)mapping);
572:   PetscOptionsFList("-islocaltoglobalmapping_type","ISLocalToGlobalMapping method","ISLocalToGlobalMappingSetType",ISLocalToGlobalMappingList,(char*)(((PetscObject)mapping)->type_name) ? ((PetscObject)mapping)->type_name : defaulttype,type,256,&flg);
573:   if (flg) {
574:     ISLocalToGlobalMappingSetType(mapping,type);
575:   }
576:   PetscOptionsEnd();
577:   return(0);
578: }

580: /*@
581:    ISLocalToGlobalMappingDestroy - Destroys a mapping between a local (0 to n)
582:    ordering and a global parallel ordering.

584:    Note Collective

586:    Input Parameters:
587: .  mapping - mapping data structure

589:    Level: advanced

591: .seealso: ISLocalToGlobalMappingCreate()
592: @*/
593: PetscErrorCode  ISLocalToGlobalMappingDestroy(ISLocalToGlobalMapping *mapping)
594: {

598:   if (!*mapping) return(0);
600:   if (--((PetscObject)(*mapping))->refct > 0) {*mapping = 0;return(0);}
601:   PetscFree((*mapping)->indices);
602:   PetscFree((*mapping)->info_procs);
603:   PetscFree((*mapping)->info_numprocs);
604:   if ((*mapping)->info_indices) {
605:     PetscInt i;

607:     PetscFree(((*mapping)->info_indices)[0]);
608:     for (i=1; i<(*mapping)->info_nproc; i++) {
609:       PetscFree(((*mapping)->info_indices)[i]);
610:     }
611:     PetscFree((*mapping)->info_indices);
612:   }
613:   if ((*mapping)->info_nodei) {
614:     PetscFree(((*mapping)->info_nodei)[0]);
615:   }
616:   PetscFree2((*mapping)->info_nodec,(*mapping)->info_nodei);
617:   if ((*mapping)->ops->destroy) {
618:     (*(*mapping)->ops->destroy)(*mapping);
619:   }
620:   PetscHeaderDestroy(mapping);
621:   *mapping = 0;
622:   return(0);
623: }

625: /*@
626:     ISLocalToGlobalMappingApplyIS - Creates from an IS in the local numbering
627:     a new index set using the global numbering defined in an ISLocalToGlobalMapping
628:     context.

630:     Collective on is

632:     Input Parameters:
633: +   mapping - mapping between local and global numbering
634: -   is - index set in local numbering

636:     Output Parameters:
637: .   newis - index set in global numbering

639:     Notes:
640:     The output IS will have the same communicator of the input IS.

642:     Level: advanced

644: .seealso: ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingCreate(),
645:           ISLocalToGlobalMappingDestroy(), ISGlobalToLocalMappingApply()
646: @*/
647: PetscErrorCode  ISLocalToGlobalMappingApplyIS(ISLocalToGlobalMapping mapping,IS is,IS *newis)
648: {
650:   PetscInt       n,*idxout;
651:   const PetscInt *idxin;


658:   ISGetLocalSize(is,&n);
659:   ISGetIndices(is,&idxin);
660:   PetscMalloc1(n,&idxout);
661:   ISLocalToGlobalMappingApply(mapping,n,idxin,idxout);
662:   ISRestoreIndices(is,&idxin);
663:   ISCreateGeneral(PetscObjectComm((PetscObject)is),n,idxout,PETSC_OWN_POINTER,newis);
664:   return(0);
665: }

667: /*@
668:    ISLocalToGlobalMappingApply - Takes a list of integers in a local numbering
669:    and converts them to the global numbering.

671:    Not collective

673:    Input Parameters:
674: +  mapping - the local to global mapping context
675: .  N - number of integers
676: -  in - input indices in local numbering

678:    Output Parameter:
679: .  out - indices in global numbering

681:    Notes:
682:    The in and out array parameters may be identical.

684:    Level: advanced

686: .seealso: ISLocalToGlobalMappingApplyBlock(), ISLocalToGlobalMappingCreate(),ISLocalToGlobalMappingDestroy(),
687:           ISLocalToGlobalMappingApplyIS(),AOCreateBasic(),AOApplicationToPetsc(),
688:           AOPetscToApplication(), ISGlobalToLocalMappingApply()

690: @*/
691: PetscErrorCode ISLocalToGlobalMappingApply(ISLocalToGlobalMapping mapping,PetscInt N,const PetscInt in[],PetscInt out[])
692: {
693:   PetscInt i,bs,Nmax;

697:   bs   = mapping->bs;
698:   Nmax = bs*mapping->n;
699:   if (bs == 1) {
700:     const PetscInt *idx = mapping->indices;
701:     for (i=0; i<N; i++) {
702:       if (in[i] < 0) {
703:         out[i] = in[i];
704:         continue;
705:       }
706:       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);
707:       out[i] = idx[in[i]];
708:     }
709:   } else {
710:     const PetscInt *idx = mapping->indices;
711:     for (i=0; i<N; i++) {
712:       if (in[i] < 0) {
713:         out[i] = in[i];
714:         continue;
715:       }
716:       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);
717:       out[i] = idx[in[i]/bs]*bs + (in[i] % bs);
718:     }
719:   }
720:   return(0);
721: }

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

726:    Not collective

728:    Input Parameters:
729: +  mapping - the local to global mapping context
730: .  N - number of integers
731: -  in - input indices in local block numbering

733:    Output Parameter:
734: .  out - indices in global block numbering

736:    Notes:
737:    The in and out array parameters may be identical.

739:    Example:
740:      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
741:      (the first block) would produce 0 and the mapping applied to 1 (the second block) would produce 3.

743:    Level: advanced

745: .seealso: ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingCreate(),ISLocalToGlobalMappingDestroy(),
746:           ISLocalToGlobalMappingApplyIS(),AOCreateBasic(),AOApplicationToPetsc(),
747:           AOPetscToApplication(), ISGlobalToLocalMappingApply()

749: @*/
750: PetscErrorCode ISLocalToGlobalMappingApplyBlock(ISLocalToGlobalMapping mapping,PetscInt N,const PetscInt in[],PetscInt out[])
751: {

755:   {
756:     PetscInt i,Nmax = mapping->n;
757:     const PetscInt *idx = mapping->indices;

759:     for (i=0; i<N; i++) {
760:       if (in[i] < 0) {
761:         out[i] = in[i];
762:         continue;
763:       }
764:       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);
765:       out[i] = idx[in[i]];
766:     }
767:   }
768:   return(0);
769: }

771: /*@
772:     ISGlobalToLocalMappingApply - Provides the local numbering for a list of integers
773:     specified with a global numbering.

775:     Not collective

777:     Input Parameters:
778: +   mapping - mapping between local and global numbering
779: .   type - IS_GTOLM_MASK - maps global indices with no local value to -1 in the output list (i.e., mask them)
780:            IS_GTOLM_DROP - drops the indices with no local value from the output list
781: .   n - number of global indices to map
782: -   idx - global indices to map

784:     Output Parameters:
785: +   nout - number of indices in output array (if type == IS_GTOLM_MASK then nout = n)
786: -   idxout - local index of each global index, one must pass in an array long enough
787:              to hold all the indices. You can call ISGlobalToLocalMappingApply() with
788:              idxout == NULL to determine the required length (returned in nout)
789:              and then allocate the required space and call ISGlobalToLocalMappingApply()
790:              a second time to set the values.

792:     Notes:
793:     Either nout or idxout may be NULL. idx and idxout may be identical.

795:     For "small" problems when using ISGlobalToLocalMappingApply() and ISGlobalToLocalMappingApplyBlock(), the ISLocalToGlobalMappingType of ISLOCALTOGLOBALMAPPINGBASIC will be used;
796:     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.
797:     Use ISLocalToGlobalMappingSetType() or call ISLocalToGlobalMappingSetFromOptions() with the option -islocaltoglobalmapping_type <basic,hash> to control which is used.

799:     Level: advanced

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

804: .seealso: ISLocalToGlobalMappingApply(), ISGlobalToLocalMappingApplyBlock(), ISLocalToGlobalMappingCreate(),
805:           ISLocalToGlobalMappingDestroy()
806: @*/
807: PetscErrorCode  ISGlobalToLocalMappingApply(ISLocalToGlobalMapping mapping,ISGlobalToLocalMappingMode type,PetscInt n,const PetscInt idx[],PetscInt *nout,PetscInt idxout[])
808: {

813:   if (!mapping->data) {
814:     ISGlobalToLocalMappingSetUp(mapping);
815:   }
816:   (*mapping->ops->globaltolocalmappingapply)(mapping,type,n,idx,nout,idxout);
817:   return(0);
818: }

820: /*@
821:     ISGlobalToLocalMappingApplyIS - Creates from an IS in the global numbering
822:     a new index set using the local numbering defined in an ISLocalToGlobalMapping
823:     context.

825:     Not collective

827:     Input Parameters:
828: +   mapping - mapping between local and global numbering
829: .   type - IS_GTOLM_MASK - maps global indices with no local value to -1 in the output list (i.e., mask them)
830:            IS_GTOLM_DROP - drops the indices with no local value from the output list
831: -   is - index set in global numbering

833:     Output Parameters:
834: .   newis - index set in local numbering

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

839:     Level: advanced

841: .seealso: ISGlobalToLocalMappingApply(), ISLocalToGlobalMappingCreate(),
842:           ISLocalToGlobalMappingDestroy()
843: @*/
844: PetscErrorCode  ISGlobalToLocalMappingApplyIS(ISLocalToGlobalMapping mapping,ISGlobalToLocalMappingMode type,IS is,IS *newis)
845: {
847:   PetscInt       n,nout,*idxout;
848:   const PetscInt *idxin;


855:   ISGetLocalSize(is,&n);
856:   ISGetIndices(is,&idxin);
857:   if (type == IS_GTOLM_MASK) {
858:     PetscMalloc1(n,&idxout);
859:   } else {
860:     ISGlobalToLocalMappingApply(mapping,type,n,idxin,&nout,NULL);
861:     PetscMalloc1(nout,&idxout);
862:   }
863:   ISGlobalToLocalMappingApply(mapping,type,n,idxin,&nout,idxout);
864:   ISRestoreIndices(is,&idxin);
865:   ISCreateGeneral(PETSC_COMM_SELF,nout,idxout,PETSC_OWN_POINTER,newis);
866:   return(0);
867: }

869: /*@
870:     ISGlobalToLocalMappingApplyBlock - Provides the local block numbering for a list of integers
871:     specified with a block global numbering.

873:     Not collective

875:     Input Parameters:
876: +   mapping - mapping between local and global numbering
877: .   type - IS_GTOLM_MASK - maps global indices with no local value to -1 in the output list (i.e., mask them)
878:            IS_GTOLM_DROP - drops the indices with no local value from the output list
879: .   n - number of global indices to map
880: -   idx - global indices to map

882:     Output Parameters:
883: +   nout - number of indices in output array (if type == IS_GTOLM_MASK then nout = n)
884: -   idxout - local index of each global index, one must pass in an array long enough
885:              to hold all the indices. You can call ISGlobalToLocalMappingApplyBlock() with
886:              idxout == NULL to determine the required length (returned in nout)
887:              and then allocate the required space and call ISGlobalToLocalMappingApplyBlock()
888:              a second time to set the values.

890:     Notes:
891:     Either nout or idxout may be NULL. idx and idxout may be identical.

893:     For "small" problems when using ISGlobalToLocalMappingApply() and ISGlobalToLocalMappingApplyBlock(), the ISLocalToGlobalMappingType of ISLOCALTOGLOBALMAPPINGBASIC will be used;
894:     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.
895:     Use ISLocalToGlobalMappingSetType() or call ISLocalToGlobalMappingSetFromOptions() with the option -islocaltoglobalmapping_type <basic,hash> to control which is used.

897:     Level: advanced

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

902: .seealso: ISLocalToGlobalMappingApply(), ISGlobalToLocalMappingApply(), ISLocalToGlobalMappingCreate(),
903:           ISLocalToGlobalMappingDestroy()
904: @*/
905: PetscErrorCode  ISGlobalToLocalMappingApplyBlock(ISLocalToGlobalMapping mapping,ISGlobalToLocalMappingMode type,
906:                                                  PetscInt n,const PetscInt idx[],PetscInt *nout,PetscInt idxout[])
907: {

912:   if (!mapping->data) {
913:     ISGlobalToLocalMappingSetUp(mapping);
914:   }
915:   (*mapping->ops->globaltolocalmappingapplyblock)(mapping,type,n,idx,nout,idxout);
916:   return(0);
917: }


920: /*@C
921:     ISLocalToGlobalMappingGetBlockInfo - Gets the neighbor information for each processor and
922:      each index shared by more than one processor

924:     Collective on ISLocalToGlobalMapping

926:     Input Parameters:
927: .   mapping - the mapping from local to global indexing

929:     Output Parameter:
930: +   nproc - number of processors that are connected to this one
931: .   proc - neighboring processors
932: .   numproc - number of indices for each subdomain (processor)
933: -   indices - indices of nodes (in local numbering) shared with neighbors (sorted by global numbering)

935:     Level: advanced

937:     Fortran Usage:
938: $        ISLocalToGlobalMpngGetInfoSize(ISLocalToGlobalMapping,PetscInt nproc,PetscInt numprocmax,ierr) followed by
939: $        ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping,PetscInt nproc, PetscInt procs[nproc],PetscInt numprocs[nproc],
940:           PetscInt indices[nproc][numprocmax],ierr)
941:         There is no ISLocalToGlobalMappingRestoreInfo() in Fortran. You must make sure that procs[], numprocs[] and
942:         indices[][] are large enough arrays, either by allocating them dynamically or defining static ones large enough.


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

954:   if (mapping->info_cached) {
955:     *nproc    = mapping->info_nproc;
956:     *procs    = mapping->info_procs;
957:     *numprocs = mapping->info_numprocs;
958:     *indices  = mapping->info_indices;
959:   } else {
960:     ISLocalToGlobalMappingGetBlockInfo_Private(mapping,nproc,procs,numprocs,indices);
961:   }
962:   return(0);
963: }

965: static PetscErrorCode  ISLocalToGlobalMappingGetBlockInfo_Private(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
966: {
968:   PetscMPIInt    size,rank,tag1,tag2,tag3,*len,*source,imdex;
969:   PetscInt       i,n = mapping->n,Ng,ng,max = 0,*lindices = mapping->indices;
970:   PetscInt       *nprocs,*owner,nsends,*sends,j,*starts,nmax,nrecvs,*recvs,proc;
971:   PetscInt       cnt,scale,*ownedsenders,*nownedsenders,rstart,nowned;
972:   PetscInt       node,nownedm,nt,*sends2,nsends2,*starts2,*lens2,*dest,nrecvs2,*starts3,*recvs2,k,*bprocs,*tmp;
973:   PetscInt       first_procs,first_numprocs,*first_indices;
974:   MPI_Request    *recv_waits,*send_waits;
975:   MPI_Status     recv_status,*send_status,*recv_statuses;
976:   MPI_Comm       comm;
977:   PetscBool      debug = PETSC_FALSE;

980:   PetscObjectGetComm((PetscObject)mapping,&comm);
981:   MPI_Comm_size(comm,&size);
982:   MPI_Comm_rank(comm,&rank);
983:   if (size == 1) {
984:     *nproc         = 0;
985:     *procs         = NULL;
986:     PetscNew(numprocs);
987:     (*numprocs)[0] = 0;
988:     PetscNew(indices);
989:     (*indices)[0]  = NULL;
990:     /* save info for reuse */
991:     mapping->info_nproc = *nproc;
992:     mapping->info_procs = *procs;
993:     mapping->info_numprocs = *numprocs;
994:     mapping->info_indices = *indices;
995:     mapping->info_cached = PETSC_TRUE;
996:     return(0);
997:   }

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

1001:   /*
1002:     Notes on ISLocalToGlobalMappingGetBlockInfo

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

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

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

1012:     nontrivial locally owned node - node that is not in the interior (i.e. has more than one
1013:            local subdomain
1014:   */
1015:   PetscObjectGetNewTag((PetscObject)mapping,&tag1);
1016:   PetscObjectGetNewTag((PetscObject)mapping,&tag2);
1017:   PetscObjectGetNewTag((PetscObject)mapping,&tag3);

1019:   for (i=0; i<n; i++) {
1020:     if (lindices[i] > max) max = lindices[i];
1021:   }
1022:   MPIU_Allreduce(&max,&Ng,1,MPIU_INT,MPI_MAX,comm);
1023:   Ng++;
1024:   MPI_Comm_size(comm,&size);
1025:   MPI_Comm_rank(comm,&rank);
1026:   scale  = Ng/size + 1;
1027:   ng     = scale; if (rank == size-1) ng = Ng - scale*(size-1); ng = PetscMax(1,ng);
1028:   rstart = scale*rank;

1030:   /* determine ownership ranges of global indices */
1031:   PetscMalloc1(2*size,&nprocs);
1032:   PetscArrayzero(nprocs,2*size);

1034:   /* determine owners of each local node  */
1035:   PetscMalloc1(n,&owner);
1036:   for (i=0; i<n; i++) {
1037:     proc             = lindices[i]/scale; /* processor that globally owns this index */
1038:     nprocs[2*proc+1] = 1;                 /* processor globally owns at least one of ours */
1039:     owner[i]         = proc;
1040:     nprocs[2*proc]++;                     /* count of how many that processor globally owns of ours */
1041:   }
1042:   nsends = 0; for (i=0; i<size; i++) nsends += nprocs[2*i+1];
1043:   PetscInfo1(mapping,"Number of global owners for my local data %D\n",nsends);

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

1049:   /* post receives for owned rows */
1050:   PetscMalloc1((2*nrecvs+1)*(nmax+1),&recvs);
1051:   PetscMalloc1(nrecvs+1,&recv_waits);
1052:   for (i=0; i<nrecvs; i++) {
1053:     MPI_Irecv(recvs+2*nmax*i,2*nmax,MPIU_INT,MPI_ANY_SOURCE,tag1,comm,recv_waits+i);
1054:   }

1056:   /* pack messages containing lists of local nodes to owners */
1057:   PetscMalloc1(2*n+1,&sends);
1058:   PetscMalloc1(size+1,&starts);
1059:   starts[0] = 0;
1060:   for (i=1; i<size; i++) starts[i] = starts[i-1] + 2*nprocs[2*i-2];
1061:   for (i=0; i<n; i++) {
1062:     sends[starts[owner[i]]++] = lindices[i];
1063:     sends[starts[owner[i]]++] = i;
1064:   }
1065:   PetscFree(owner);
1066:   starts[0] = 0;
1067:   for (i=1; i<size; i++) starts[i] = starts[i-1] + 2*nprocs[2*i-2];

1069:   /* send the messages */
1070:   PetscMalloc1(nsends+1,&send_waits);
1071:   PetscMalloc1(nsends+1,&dest);
1072:   cnt = 0;
1073:   for (i=0; i<size; i++) {
1074:     if (nprocs[2*i]) {
1075:       MPI_Isend(sends+starts[i],2*nprocs[2*i],MPIU_INT,i,tag1,comm,send_waits+cnt);
1076:       dest[cnt] = i;
1077:       cnt++;
1078:     }
1079:   }
1080:   PetscFree(starts);

1082:   /* wait on receives */
1083:   PetscMalloc1(nrecvs+1,&source);
1084:   PetscMalloc1(nrecvs+1,&len);
1085:   cnt  = nrecvs;
1086:   PetscCalloc1(ng+1,&nownedsenders);
1087:   while (cnt) {
1088:     MPI_Waitany(nrecvs,recv_waits,&imdex,&recv_status);
1089:     /* unpack receives into our local space */
1090:     MPI_Get_count(&recv_status,MPIU_INT,&len[imdex]);
1091:     source[imdex] = recv_status.MPI_SOURCE;
1092:     len[imdex]    = len[imdex]/2;
1093:     /* count how many local owners for each of my global owned indices */
1094:     for (i=0; i<len[imdex]; i++) nownedsenders[recvs[2*imdex*nmax+2*i]-rstart]++;
1095:     cnt--;
1096:   }
1097:   PetscFree(recv_waits);

1099:   /* count how many globally owned indices are on an edge multiplied by how many processors own them. */
1100:   nowned  = 0;
1101:   nownedm = 0;
1102:   for (i=0; i<ng; i++) {
1103:     if (nownedsenders[i] > 1) {nownedm += nownedsenders[i]; nowned++;}
1104:   }

1106:   /* create single array to contain rank of all local owners of each globally owned index */
1107:   PetscMalloc1(nownedm+1,&ownedsenders);
1108:   PetscMalloc1(ng+1,&starts);
1109:   starts[0] = 0;
1110:   for (i=1; i<ng; i++) {
1111:     if (nownedsenders[i-1] > 1) starts[i] = starts[i-1] + nownedsenders[i-1];
1112:     else starts[i] = starts[i-1];
1113:   }

1115:   /* for each nontrival globally owned node list all arriving processors */
1116:   for (i=0; i<nrecvs; i++) {
1117:     for (j=0; j<len[i]; j++) {
1118:       node = recvs[2*i*nmax+2*j]-rstart;
1119:       if (nownedsenders[node] > 1) ownedsenders[starts[node]++] = source[i];
1120:     }
1121:   }

1123:   if (debug) { /* -----------------------------------  */
1124:     starts[0] = 0;
1125:     for (i=1; i<ng; i++) {
1126:       if (nownedsenders[i-1] > 1) starts[i] = starts[i-1] + nownedsenders[i-1];
1127:       else starts[i] = starts[i-1];
1128:     }
1129:     for (i=0; i<ng; i++) {
1130:       if (nownedsenders[i] > 1) {
1131:         PetscSynchronizedPrintf(comm,"[%d] global node %D local owner processors: ",rank,i+rstart);
1132:         for (j=0; j<nownedsenders[i]; j++) {
1133:           PetscSynchronizedPrintf(comm,"%D ",ownedsenders[starts[i]+j]);
1134:         }
1135:         PetscSynchronizedPrintf(comm,"\n");
1136:       }
1137:     }
1138:     PetscSynchronizedFlush(comm,PETSC_STDOUT);
1139:   } /* -----------------------------------  */

1141:   /* wait on original sends */
1142:   if (nsends) {
1143:     PetscMalloc1(nsends,&send_status);
1144:     MPI_Waitall(nsends,send_waits,send_status);
1145:     PetscFree(send_status);
1146:   }
1147:   PetscFree(send_waits);
1148:   PetscFree(sends);
1149:   PetscFree(nprocs);

1151:   /* pack messages to send back to local owners */
1152:   starts[0] = 0;
1153:   for (i=1; i<ng; i++) {
1154:     if (nownedsenders[i-1] > 1) starts[i] = starts[i-1] + nownedsenders[i-1];
1155:     else starts[i] = starts[i-1];
1156:   }
1157:   nsends2 = nrecvs;
1158:   PetscMalloc1(nsends2+1,&nprocs); /* length of each message */
1159:   for (i=0; i<nrecvs; i++) {
1160:     nprocs[i] = 1;
1161:     for (j=0; j<len[i]; j++) {
1162:       node = recvs[2*i*nmax+2*j]-rstart;
1163:       if (nownedsenders[node] > 1) nprocs[i] += 2 + nownedsenders[node];
1164:     }
1165:   }
1166:   nt = 0;
1167:   for (i=0; i<nsends2; i++) nt += nprocs[i];

1169:   PetscMalloc1(nt+1,&sends2);
1170:   PetscMalloc1(nsends2+1,&starts2);

1172:   starts2[0] = 0;
1173:   for (i=1; i<nsends2; i++) starts2[i] = starts2[i-1] + nprocs[i-1];
1174:   /*
1175:      Each message is 1 + nprocs[i] long, and consists of
1176:        (0) the number of nodes being sent back
1177:        (1) the local node number,
1178:        (2) the number of processors sharing it,
1179:        (3) the processors sharing it
1180:   */
1181:   for (i=0; i<nsends2; i++) {
1182:     cnt = 1;
1183:     sends2[starts2[i]] = 0;
1184:     for (j=0; j<len[i]; j++) {
1185:       node = recvs[2*i*nmax+2*j]-rstart;
1186:       if (nownedsenders[node] > 1) {
1187:         sends2[starts2[i]]++;
1188:         sends2[starts2[i]+cnt++] = recvs[2*i*nmax+2*j+1];
1189:         sends2[starts2[i]+cnt++] = nownedsenders[node];
1190:         PetscArraycpy(&sends2[starts2[i]+cnt],&ownedsenders[starts[node]],nownedsenders[node]);
1191:         cnt += nownedsenders[node];
1192:       }
1193:     }
1194:   }

1196:   /* receive the message lengths */
1197:   nrecvs2 = nsends;
1198:   PetscMalloc1(nrecvs2+1,&lens2);
1199:   PetscMalloc1(nrecvs2+1,&starts3);
1200:   PetscMalloc1(nrecvs2+1,&recv_waits);
1201:   for (i=0; i<nrecvs2; i++) {
1202:     MPI_Irecv(&lens2[i],1,MPIU_INT,dest[i],tag2,comm,recv_waits+i);
1203:   }

1205:   /* send the message lengths */
1206:   for (i=0; i<nsends2; i++) {
1207:     MPI_Send(&nprocs[i],1,MPIU_INT,source[i],tag2,comm);
1208:   }

1210:   /* wait on receives of lens */
1211:   if (nrecvs2) {
1212:     PetscMalloc1(nrecvs2,&recv_statuses);
1213:     MPI_Waitall(nrecvs2,recv_waits,recv_statuses);
1214:     PetscFree(recv_statuses);
1215:   }
1216:   PetscFree(recv_waits);

1218:   starts3[0] = 0;
1219:   nt         = 0;
1220:   for (i=0; i<nrecvs2-1; i++) {
1221:     starts3[i+1] = starts3[i] + lens2[i];
1222:     nt          += lens2[i];
1223:   }
1224:   if (nrecvs2) nt += lens2[nrecvs2-1];

1226:   PetscMalloc1(nt+1,&recvs2);
1227:   PetscMalloc1(nrecvs2+1,&recv_waits);
1228:   for (i=0; i<nrecvs2; i++) {
1229:     MPI_Irecv(recvs2+starts3[i],lens2[i],MPIU_INT,dest[i],tag3,comm,recv_waits+i);
1230:   }

1232:   /* send the messages */
1233:   PetscMalloc1(nsends2+1,&send_waits);
1234:   for (i=0; i<nsends2; i++) {
1235:     MPI_Isend(sends2+starts2[i],nprocs[i],MPIU_INT,source[i],tag3,comm,send_waits+i);
1236:   }

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

1247:   if (debug) { /* -----------------------------------  */
1248:     cnt = 0;
1249:     for (i=0; i<nrecvs2; i++) {
1250:       nt = recvs2[cnt++];
1251:       for (j=0; j<nt; j++) {
1252:         PetscSynchronizedPrintf(comm,"[%d] local node %D number of subdomains %D: ",rank,recvs2[cnt],recvs2[cnt+1]);
1253:         for (k=0; k<recvs2[cnt+1]; k++) {
1254:           PetscSynchronizedPrintf(comm,"%D ",recvs2[cnt+2+k]);
1255:         }
1256:         cnt += 2 + recvs2[cnt+1];
1257:         PetscSynchronizedPrintf(comm,"\n");
1258:       }
1259:     }
1260:     PetscSynchronizedFlush(comm,PETSC_STDOUT);
1261:   } /* -----------------------------------  */

1263:   /* count number subdomains for each local node */
1264:   PetscCalloc1(size,&nprocs);
1265:   cnt  = 0;
1266:   for (i=0; i<nrecvs2; i++) {
1267:     nt = recvs2[cnt++];
1268:     for (j=0; j<nt; j++) {
1269:       for (k=0; k<recvs2[cnt+1]; k++) nprocs[recvs2[cnt+2+k]]++;
1270:       cnt += 2 + recvs2[cnt+1];
1271:     }
1272:   }
1273:   nt = 0; for (i=0; i<size; i++) nt += (nprocs[i] > 0);
1274:   *nproc    = nt;
1275:   PetscMalloc1(nt+1,procs);
1276:   PetscMalloc1(nt+1,numprocs);
1277:   PetscMalloc1(nt+1,indices);
1278:   for (i=0;i<nt+1;i++) (*indices)[i]=NULL;
1279:   PetscMalloc1(size,&bprocs);
1280:   cnt  = 0;
1281:   for (i=0; i<size; i++) {
1282:     if (nprocs[i] > 0) {
1283:       bprocs[i]        = cnt;
1284:       (*procs)[cnt]    = i;
1285:       (*numprocs)[cnt] = nprocs[i];
1286:       PetscMalloc1(nprocs[i],&(*indices)[cnt]);
1287:       cnt++;
1288:     }
1289:   }

1291:   /* make the list of subdomains for each nontrivial local node */
1292:   PetscArrayzero(*numprocs,nt);
1293:   cnt  = 0;
1294:   for (i=0; i<nrecvs2; i++) {
1295:     nt = recvs2[cnt++];
1296:     for (j=0; j<nt; j++) {
1297:       for (k=0; k<recvs2[cnt+1]; k++) (*indices)[bprocs[recvs2[cnt+2+k]]][(*numprocs)[bprocs[recvs2[cnt+2+k]]]++] = recvs2[cnt];
1298:       cnt += 2 + recvs2[cnt+1];
1299:     }
1300:   }
1301:   PetscFree(bprocs);
1302:   PetscFree(recvs2);

1304:   /* sort the node indexing by their global numbers */
1305:   nt = *nproc;
1306:   for (i=0; i<nt; i++) {
1307:     PetscMalloc1((*numprocs)[i],&tmp);
1308:     for (j=0; j<(*numprocs)[i]; j++) tmp[j] = lindices[(*indices)[i][j]];
1309:     PetscSortIntWithArray((*numprocs)[i],tmp,(*indices)[i]);
1310:     PetscFree(tmp);
1311:   }

1313:   if (debug) { /* -----------------------------------  */
1314:     nt = *nproc;
1315:     for (i=0; i<nt; i++) {
1316:       PetscSynchronizedPrintf(comm,"[%d] subdomain %D number of indices %D: ",rank,(*procs)[i],(*numprocs)[i]);
1317:       for (j=0; j<(*numprocs)[i]; j++) {
1318:         PetscSynchronizedPrintf(comm,"%D ",(*indices)[i][j]);
1319:       }
1320:       PetscSynchronizedPrintf(comm,"\n");
1321:     }
1322:     PetscSynchronizedFlush(comm,PETSC_STDOUT);
1323:   } /* -----------------------------------  */

1325:   /* wait on sends */
1326:   if (nsends2) {
1327:     PetscMalloc1(nsends2,&send_status);
1328:     MPI_Waitall(nsends2,send_waits,send_status);
1329:     PetscFree(send_status);
1330:   }

1332:   PetscFree(starts3);
1333:   PetscFree(dest);
1334:   PetscFree(send_waits);

1336:   PetscFree(nownedsenders);
1337:   PetscFree(ownedsenders);
1338:   PetscFree(starts);
1339:   PetscFree(starts2);
1340:   PetscFree(lens2);

1342:   PetscFree(source);
1343:   PetscFree(len);
1344:   PetscFree(recvs);
1345:   PetscFree(nprocs);
1346:   PetscFree(sends2);

1348:   /* put the information about myself as the first entry in the list */
1349:   first_procs    = (*procs)[0];
1350:   first_numprocs = (*numprocs)[0];
1351:   first_indices  = (*indices)[0];
1352:   for (i=0; i<*nproc; i++) {
1353:     if ((*procs)[i] == rank) {
1354:       (*procs)[0]    = (*procs)[i];
1355:       (*numprocs)[0] = (*numprocs)[i];
1356:       (*indices)[0]  = (*indices)[i];
1357:       (*procs)[i]    = first_procs;
1358:       (*numprocs)[i] = first_numprocs;
1359:       (*indices)[i]  = first_indices;
1360:       break;
1361:     }
1362:   }

1364:   /* save info for reuse */
1365:   mapping->info_nproc = *nproc;
1366:   mapping->info_procs = *procs;
1367:   mapping->info_numprocs = *numprocs;
1368:   mapping->info_indices = *indices;
1369:   mapping->info_cached = PETSC_TRUE;
1370:   return(0);
1371: }

1373: /*@C
1374:     ISLocalToGlobalMappingRestoreBlockInfo - Frees the memory allocated by ISLocalToGlobalMappingGetBlockInfo()

1376:     Collective on ISLocalToGlobalMapping

1378:     Input Parameters:
1379: .   mapping - the mapping from local to global indexing

1381:     Output Parameter:
1382: +   nproc - number of processors that are connected to this one
1383: .   proc - neighboring processors
1384: .   numproc - number of indices for each processor
1385: -   indices - indices of local nodes shared with neighbor (sorted by global numbering)

1387:     Level: advanced

1389: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
1390:           ISLocalToGlobalMappingGetInfo()
1391: @*/
1392: PetscErrorCode  ISLocalToGlobalMappingRestoreBlockInfo(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
1393: {

1398:   if (mapping->info_free) {
1399:     PetscFree(*numprocs);
1400:     if (*indices) {
1401:       PetscInt i;

1403:       PetscFree((*indices)[0]);
1404:       for (i=1; i<*nproc; i++) {
1405:         PetscFree((*indices)[i]);
1406:       }
1407:       PetscFree(*indices);
1408:     }
1409:   }
1410:   *nproc    = 0;
1411:   *procs    = NULL;
1412:   *numprocs = NULL;
1413:   *indices  = NULL;
1414:   return(0);
1415: }

1417: /*@C
1418:     ISLocalToGlobalMappingGetInfo - Gets the neighbor information for each processor and
1419:      each index shared by more than one processor

1421:     Collective on ISLocalToGlobalMapping

1423:     Input Parameters:
1424: .   mapping - the mapping from local to global indexing

1426:     Output Parameter:
1427: +   nproc - number of processors that are connected to this one
1428: .   proc - neighboring processors
1429: .   numproc - number of indices for each subdomain (processor)
1430: -   indices - indices of nodes (in local numbering) shared with neighbors (sorted by global numbering)

1432:     Level: advanced

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

1436:     Fortran Usage:
1437: $        ISLocalToGlobalMpngGetInfoSize(ISLocalToGlobalMapping,PetscInt nproc,PetscInt numprocmax,ierr) followed by
1438: $        ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping,PetscInt nproc, PetscInt procs[nproc],PetscInt numprocs[nproc],
1439:           PetscInt indices[nproc][numprocmax],ierr)
1440:         There is no ISLocalToGlobalMappingRestoreInfo() in Fortran. You must make sure that procs[], numprocs[] and
1441:         indices[][] are large enough arrays, either by allocating them dynamically or defining static ones large enough.


1444: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
1445:           ISLocalToGlobalMappingRestoreInfo()
1446: @*/
1447: PetscErrorCode  ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
1448: {
1450:   PetscInt       **bindices = NULL,*bnumprocs = NULL,bs = mapping->bs,i,j,k;

1454:   ISLocalToGlobalMappingGetBlockInfo(mapping,nproc,procs,&bnumprocs,&bindices);
1455:   if (bs > 1) { /* we need to expand the cached info */
1456:     PetscCalloc1(*nproc,&*indices);
1457:     PetscCalloc1(*nproc,&*numprocs);
1458:     for (i=0; i<*nproc; i++) {
1459:       PetscMalloc1(bs*bnumprocs[i],&(*indices)[i]);
1460:       for (j=0; j<bnumprocs[i]; j++) {
1461:         for (k=0; k<bs; k++) {
1462:           (*indices)[i][j*bs+k] = bs*bindices[i][j] + k;
1463:         }
1464:       }
1465:       (*numprocs)[i] = bnumprocs[i]*bs;
1466:     }
1467:     mapping->info_free = PETSC_TRUE;
1468:   } else {
1469:     *numprocs = bnumprocs;
1470:     *indices  = bindices;
1471:   }
1472:   return(0);
1473: }

1475: /*@C
1476:     ISLocalToGlobalMappingRestoreInfo - Frees the memory allocated by ISLocalToGlobalMappingGetInfo()

1478:     Collective on ISLocalToGlobalMapping

1480:     Input Parameters:
1481: .   mapping - the mapping from local to global indexing

1483:     Output Parameter:
1484: +   nproc - number of processors that are connected to this one
1485: .   proc - neighboring processors
1486: .   numproc - number of indices for each processor
1487: -   indices - indices of local nodes shared with neighbor (sorted by global numbering)

1489:     Level: advanced

1491: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
1492:           ISLocalToGlobalMappingGetInfo()
1493: @*/
1494: PetscErrorCode  ISLocalToGlobalMappingRestoreInfo(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
1495: {

1499:   ISLocalToGlobalMappingRestoreBlockInfo(mapping,nproc,procs,numprocs,indices);
1500:   return(0);
1501: }

1503: /*@C
1504:     ISLocalToGlobalMappingGetNodeInfo - Gets the neighbor information for each node

1506:     Collective on ISLocalToGlobalMapping

1508:     Input Parameters:
1509: .   mapping - the mapping from local to global indexing

1511:     Output Parameter:
1512: +   nnodes - number of local nodes (same ISLocalToGlobalMappingGetSize())
1513: .   count - number of neighboring processors per node
1514: -   indices - indices of processes sharing the node (sorted)

1516:     Level: advanced

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

1520: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
1521:           ISLocalToGlobalMappingGetInfo(), ISLocalToGlobalMappingRestoreNodeInfo()
1522: @*/
1523: PetscErrorCode  ISLocalToGlobalMappingGetNodeInfo(ISLocalToGlobalMapping mapping,PetscInt *nnodes,PetscInt *count[],PetscInt **indices[])
1524: {
1525:   PetscInt       n;

1530:   ISLocalToGlobalMappingGetSize(mapping,&n);
1531:   if (!mapping->info_nodec) {
1532:     PetscInt i,m,n_neigh,*neigh,*n_shared,**shared;

1534:     PetscMalloc2(n+1,&mapping->info_nodec,n,&mapping->info_nodei);
1535:     ISLocalToGlobalMappingGetInfo(mapping,&n_neigh,&neigh,&n_shared,&shared);
1536:     for (i=0;i<n;i++) { mapping->info_nodec[i] = 1;}
1537:     m = n;
1538:     mapping->info_nodec[n] = 0;
1539:     for (i=1;i<n_neigh;i++) {
1540:       PetscInt j;

1542:       m += n_shared[i];
1543:       for (j=0;j<n_shared[i];j++) mapping->info_nodec[shared[i][j]] += 1;
1544:     }
1545:     if (n) { PetscMalloc1(m,&mapping->info_nodei[0]); }
1546:     for (i=1;i<n;i++) mapping->info_nodei[i] = mapping->info_nodei[i-1] + mapping->info_nodec[i-1];
1547:     PetscArrayzero(mapping->info_nodec,n);
1548:     for (i=0;i<n;i++) { mapping->info_nodec[i] = 1; mapping->info_nodei[i][0] = neigh[0]; }
1549:     for (i=1;i<n_neigh;i++) {
1550:       PetscInt j;

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

1555:         mapping->info_nodei[k][mapping->info_nodec[k]] = neigh[i];
1556:         mapping->info_nodec[k] += 1;
1557:       }
1558:     }
1559:     for (i=0;i<n;i++) { PetscSortRemoveDupsInt(&mapping->info_nodec[i],mapping->info_nodei[i]); }
1560:     ISLocalToGlobalMappingRestoreInfo(mapping,&n_neigh,&neigh,&n_shared,&shared);
1561:   }
1562:   if (nnodes)  *nnodes  = n;
1563:   if (count)   *count   = mapping->info_nodec;
1564:   if (indices) *indices = mapping->info_nodei;
1565:   return(0);
1566: }

1568: /*@C
1569:     ISLocalToGlobalMappingRestoreNodeInfo - Frees the memory allocated by ISLocalToGlobalMappingGetNodeInfo()

1571:     Collective on ISLocalToGlobalMapping

1573:     Input Parameters:
1574: .   mapping - the mapping from local to global indexing

1576:     Output Parameter:
1577: +   nnodes - number of local nodes
1578: .   count - number of neighboring processors per node
1579: -   indices - indices of processes sharing the node (sorted)

1581:     Level: advanced

1583: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
1584:           ISLocalToGlobalMappingGetInfo()
1585: @*/
1586: PetscErrorCode  ISLocalToGlobalMappingRestoreNodeInfo(ISLocalToGlobalMapping mapping,PetscInt *nnodes,PetscInt *count[],PetscInt **indices[])
1587: {
1590:   if (nnodes)  *nnodes  = 0;
1591:   if (count)   *count   = NULL;
1592:   if (indices) *indices = NULL;
1593:   return(0);
1594: }

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

1599:    Not Collective

1601:    Input Arguments:
1602: . ltog - local to global mapping

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

1607:    Level: advanced

1609:    Notes:
1610:     ISLocalToGlobalMappingGetSize() returns the length the this array

1612: .seealso: ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingRestoreIndices(), ISLocalToGlobalMappingGetBlockIndices(), ISLocalToGlobalMappingRestoreBlockIndices()
1613: @*/
1614: PetscErrorCode  ISLocalToGlobalMappingGetIndices(ISLocalToGlobalMapping ltog,const PetscInt **array)
1615: {
1619:   if (ltog->bs == 1) {
1620:     *array = ltog->indices;
1621:   } else {
1622:     PetscInt       *jj,k,i,j,n = ltog->n, bs = ltog->bs;
1623:     const PetscInt *ii;

1626:     PetscMalloc1(bs*n,&jj);
1627:     *array = jj;
1628:     k    = 0;
1629:     ii   = ltog->indices;
1630:     for (i=0; i<n; i++)
1631:       for (j=0; j<bs; j++)
1632:         jj[k++] = bs*ii[i] + j;
1633:   }
1634:   return(0);
1635: }

1637: /*@C
1638:    ISLocalToGlobalMappingRestoreIndices - Restore indices obtained with ISLocalToGlobalMappingGetIndices()

1640:    Not Collective

1642:    Input Arguments:
1643: + ltog - local to global mapping
1644: - array - array of indices

1646:    Level: advanced

1648: .seealso: ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingGetIndices()
1649: @*/
1650: PetscErrorCode  ISLocalToGlobalMappingRestoreIndices(ISLocalToGlobalMapping ltog,const PetscInt **array)
1651: {
1655:   if (ltog->bs == 1 && *array != ltog->indices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_BADPTR,"Trying to return mismatched pointer");

1657:   if (ltog->bs > 1) {
1659:     PetscFree(*(void**)array);
1660:   }
1661:   return(0);
1662: }

1664: /*@C
1665:    ISLocalToGlobalMappingGetBlockIndices - Get global indices for every local block

1667:    Not Collective

1669:    Input Arguments:
1670: . ltog - local to global mapping

1672:    Output Arguments:
1673: . array - array of indices

1675:    Level: advanced

1677: .seealso: ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingRestoreBlockIndices()
1678: @*/
1679: PetscErrorCode  ISLocalToGlobalMappingGetBlockIndices(ISLocalToGlobalMapping ltog,const PetscInt **array)
1680: {
1684:   *array = ltog->indices;
1685:   return(0);
1686: }

1688: /*@C
1689:    ISLocalToGlobalMappingRestoreBlockIndices - Restore indices obtained with ISLocalToGlobalMappingGetBlockIndices()

1691:    Not Collective

1693:    Input Arguments:
1694: + ltog - local to global mapping
1695: - array - array of indices

1697:    Level: advanced

1699: .seealso: ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingGetIndices()
1700: @*/
1701: PetscErrorCode  ISLocalToGlobalMappingRestoreBlockIndices(ISLocalToGlobalMapping ltog,const PetscInt **array)
1702: {
1706:   if (*array != ltog->indices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_BADPTR,"Trying to return mismatched pointer");
1707:   *array = NULL;
1708:   return(0);
1709: }

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

1714:    Not Collective

1716:    Input Arguments:
1717: + comm - communicator for the new mapping, must contain the communicator of every mapping to concatenate
1718: . n - number of mappings to concatenate
1719: - ltogs - local to global mappings

1721:    Output Arguments:
1722: . ltogcat - new mapping

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

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

1728:    Level: advanced

1730: .seealso: ISLocalToGlobalMappingCreate()
1731: @*/
1732: PetscErrorCode ISLocalToGlobalMappingConcatenate(MPI_Comm comm,PetscInt n,const ISLocalToGlobalMapping ltogs[],ISLocalToGlobalMapping *ltogcat)
1733: {
1734:   PetscInt       i,cnt,m,*idx;

1738:   if (n < 0) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have a non-negative number of mappings, given %D",n);
1742:   for (cnt=0,i=0; i<n; i++) {
1743:     ISLocalToGlobalMappingGetSize(ltogs[i],&m);
1744:     cnt += m;
1745:   }
1746:   PetscMalloc1(cnt,&idx);
1747:   for (cnt=0,i=0; i<n; i++) {
1748:     const PetscInt *subidx;
1749:     ISLocalToGlobalMappingGetSize(ltogs[i],&m);
1750:     ISLocalToGlobalMappingGetIndices(ltogs[i],&subidx);
1751:     PetscArraycpy(&idx[cnt],subidx,m);
1752:     ISLocalToGlobalMappingRestoreIndices(ltogs[i],&subidx);
1753:     cnt += m;
1754:   }
1755:   ISLocalToGlobalMappingCreate(comm,1,cnt,idx,PETSC_OWN_POINTER,ltogcat);
1756:   return(0);
1757: }

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

1763:    Options Database Keys:
1764: .   -islocaltoglobalmapping_type basic - select this method

1766:    Level: beginner

1768: .seealso:  ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingSetType(), ISLOCALTOGLOBALMAPPINGHASH
1769: M*/
1770: PETSC_EXTERN PetscErrorCode ISLocalToGlobalMappingCreate_Basic(ISLocalToGlobalMapping ltog)
1771: {
1773:   ltog->ops->globaltolocalmappingapply      = ISGlobalToLocalMappingApply_Basic;
1774:   ltog->ops->globaltolocalmappingsetup      = ISGlobalToLocalMappingSetUp_Basic;
1775:   ltog->ops->globaltolocalmappingapplyblock = ISGlobalToLocalMappingApplyBlock_Basic;
1776:   ltog->ops->destroy                        = ISLocalToGlobalMappingDestroy_Basic;
1777:   return(0);
1778: }

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

1784:    Options Database Keys:
1785: .   -islocaltoglobalmapping_type hash - select this method


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

1791:    Level: beginner

1793: .seealso:  ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingSetType(), ISLOCALTOGLOBALMAPPINGHASH
1794: M*/
1795: PETSC_EXTERN PetscErrorCode ISLocalToGlobalMappingCreate_Hash(ISLocalToGlobalMapping ltog)
1796: {
1798:   ltog->ops->globaltolocalmappingapply      = ISGlobalToLocalMappingApply_Hash;
1799:   ltog->ops->globaltolocalmappingsetup      = ISGlobalToLocalMappingSetUp_Hash;
1800:   ltog->ops->globaltolocalmappingapplyblock = ISGlobalToLocalMappingApplyBlock_Hash;
1801:   ltog->ops->destroy                        = ISLocalToGlobalMappingDestroy_Hash;
1802:   return(0);
1803: }


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

1809:    Not Collective

1811:    Input Parameters:
1812: +  sname - name of a new method
1813: -  routine_create - routine to create method context

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

1818:    Sample usage:
1819: .vb
1820:    ISLocalToGlobalMappingRegister("my_mapper",MyCreate);
1821: .ve

1823:    Then, your mapping can be chosen with the procedural interface via
1824: $     ISLocalToGlobalMappingSetType(ltog,"my_mapper")
1825:    or at runtime via the option
1826: $     -islocaltoglobalmapping_type my_mapper

1828:    Level: advanced

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

1832: @*/
1833: PetscErrorCode  ISLocalToGlobalMappingRegister(const char sname[],PetscErrorCode (*function)(ISLocalToGlobalMapping))
1834: {

1838:   ISInitializePackage();
1839:   PetscFunctionListAdd(&ISLocalToGlobalMappingList,sname,function);
1840:   return(0);
1841: }

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

1846:    Logically Collective on ISLocalToGlobalMapping

1848:    Input Parameters:
1849: +  ltog - the ISLocalToGlobalMapping object
1850: -  type - a known method

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

1856:    Notes:
1857:    See "petsc/include/petscis.h" for available methods

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

1863:   Level: intermediate

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

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

1870: @*/
1871: PetscErrorCode  ISLocalToGlobalMappingSetType(ISLocalToGlobalMapping ltog, ISLocalToGlobalMappingType type)
1872: {
1873:   PetscErrorCode ierr,(*r)(ISLocalToGlobalMapping);
1874:   PetscBool      match;


1880:   PetscObjectTypeCompare((PetscObject)ltog,type,&match);
1881:   if (match) return(0);

1883:    PetscFunctionListFind(ISLocalToGlobalMappingList,type,&r);
1884:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)ltog),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested ISLocalToGlobalMapping type %s",type);
1885:   /* Destroy the previous private LTOG context */
1886:   if (ltog->ops->destroy) {
1887:     (*ltog->ops->destroy)(ltog);
1888:     ltog->ops->destroy = NULL;
1889:   }
1890:   PetscObjectChangeTypeName((PetscObject)ltog,type);
1891:   (*r)(ltog);
1892:   return(0);
1893: }

1895: PetscBool ISLocalToGlobalMappingRegisterAllCalled = PETSC_FALSE;

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

1900:   Not Collective

1902:   Level: advanced

1904: .seealso:  ISRegister(),  ISLocalToGlobalRegister()
1905: @*/
1906: PetscErrorCode  ISLocalToGlobalMappingRegisterAll(void)
1907: {

1911:   if (ISLocalToGlobalMappingRegisterAllCalled) return(0);
1912:   ISLocalToGlobalMappingRegisterAllCalled = PETSC_TRUE;

1914:   ISLocalToGlobalMappingRegister(ISLOCALTOGLOBALMAPPINGBASIC, ISLocalToGlobalMappingCreate_Basic);
1915:   ISLocalToGlobalMappingRegister(ISLOCALTOGLOBALMAPPINGHASH, ISLocalToGlobalMappingCreate_Hash);
1916:   return(0);
1917: }