Actual source code: isltog.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  2: #include <petsc/private/isimpl.h>    /*I "petscis.h"  I*/
  3: #include <petscsf.h>
  4: #include <petscviewer.h>

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


 12: PetscErrorCode ISG2LMapApply(ISLocalToGlobalMapping mapping,PetscInt n,const PetscInt in[],PetscInt out[])
 13: {
 15:   PetscInt       i,start,end;

 18:   if (!mapping->globals) {
 19:     ISGlobalToLocalMappingApply(mapping,IS_GTOLM_MASK,0,0,0,0);
 20:   }
 21:   start = mapping->globalstart;
 22:   end = mapping->globalend;
 23:   for (i=0; i<n; i++) {
 24:     if (in[i] < 0)          out[i] = in[i];
 25:     else if (in[i] < start) out[i] = -1;
 26:     else if (in[i] > end)   out[i] = -1;
 27:     else                    out[i] = mapping->globals[in[i] - start];
 28:   }
 29:   return(0);
 30: }


 35: /*@
 36:     ISLocalToGlobalMappingGetSize - Gets the local size of a local to global mapping

 38:     Not Collective

 40:     Input Parameter:
 41: .   ltog - local to global mapping

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

 46:     Level: advanced

 48:     Concepts: mapping^local to global

 50: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate()
 51: @*/
 52: PetscErrorCode  ISLocalToGlobalMappingGetSize(ISLocalToGlobalMapping mapping,PetscInt *n)
 53: {
 57:   *n = mapping->bs*mapping->n;
 58:   return(0);
 59: }

 63: /*@C
 64:     ISLocalToGlobalMappingView - View a local to global mapping

 66:     Not Collective

 68:     Input Parameters:
 69: +   ltog - local to global mapping
 70: -   viewer - viewer

 72:     Level: advanced

 74:     Concepts: mapping^local to global

 76: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate()
 77: @*/
 78: PetscErrorCode  ISLocalToGlobalMappingView(ISLocalToGlobalMapping mapping,PetscViewer viewer)
 79: {
 80:   PetscInt       i;
 81:   PetscMPIInt    rank;
 82:   PetscBool      iascii;

 87:   if (!viewer) {
 88:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)mapping),&viewer);
 89:   }

 92:   MPI_Comm_rank(PetscObjectComm((PetscObject)mapping),&rank);
 93:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 94:   if (iascii) {
 95:     PetscObjectPrintClassNamePrefixType((PetscObject)mapping,viewer);
 96:     PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
 97:     for (i=0; i<mapping->n; i++) {
 98:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,mapping->indices[i]);
 99:     }
100:     PetscViewerFlush(viewer);
101:     PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
102:   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for ISLocalToGlobalMapping",((PetscObject)viewer)->type_name);
103:   return(0);
104: }

108: /*@
109:     ISLocalToGlobalMappingCreateIS - Creates a mapping between a local (0 to n)
110:     ordering and a global parallel ordering.

112:     Not collective

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

117:     Output Parameter:
118: .   mapping - new mapping data structure

120:     Notes: the block size of the IS determines the block size of the mapping
121:     Level: advanced

123:     Concepts: mapping^local to global

125: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate()
126: @*/
127: PetscErrorCode  ISLocalToGlobalMappingCreateIS(IS is,ISLocalToGlobalMapping *mapping)
128: {
130:   PetscInt       n,bs;
131:   const PetscInt *indices;
132:   MPI_Comm       comm;
133:   PetscBool      isblock;


139:   PetscObjectGetComm((PetscObject)is,&comm);
140:   ISGetLocalSize(is,&n);
141:   PetscObjectTypeCompare((PetscObject)is,ISBLOCK,&isblock);
142:   if (!isblock) {
143:     ISGetIndices(is,&indices);
144:     ISLocalToGlobalMappingCreate(comm,1,n,indices,PETSC_COPY_VALUES,mapping);
145:     ISRestoreIndices(is,&indices);
146:   } else {
147:     ISGetBlockSize(is,&bs);
148:     ISBlockGetIndices(is,&indices);
149:     ISLocalToGlobalMappingCreate(comm,bs,n/bs,indices,PETSC_COPY_VALUES,mapping);
150:     ISBlockRestoreIndices(is,&indices);
151:   }
152:   return(0);
153: }

157: /*@C
158:     ISLocalToGlobalMappingCreateSF - Creates a mapping between a local (0 to n)
159:     ordering and a global parallel ordering.

161:     Collective

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

167:     Output Parameter:
168: .   mapping - new mapping data structure

170:     Level: advanced

172:     Concepts: mapping^local to global

174: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingCreateIS()
175: @*/
176: PetscErrorCode ISLocalToGlobalMappingCreateSF(PetscSF sf,PetscInt start,ISLocalToGlobalMapping *mapping)
177: {
179:   PetscInt       i,maxlocal,nroots,nleaves,*globals,*ltog;
180:   const PetscInt *ilocal;
181:   MPI_Comm       comm;


187:   PetscObjectGetComm((PetscObject)sf,&comm);
188:   PetscSFGetGraph(sf,&nroots,&nleaves,&ilocal,NULL);
189:   if (ilocal) {
190:     for (i=0,maxlocal=0; i<nleaves; i++) maxlocal = PetscMax(maxlocal,ilocal[i]+1);
191:   }
192:   else maxlocal = nleaves;
193:   PetscMalloc1(nroots,&globals);
194:   PetscMalloc1(maxlocal,&ltog);
195:   for (i=0; i<nroots; i++) globals[i] = start + i;
196:   for (i=0; i<maxlocal; i++) ltog[i] = -1;
197:   PetscSFBcastBegin(sf,MPIU_INT,globals,ltog);
198:   PetscSFBcastEnd(sf,MPIU_INT,globals,ltog);
199:   ISLocalToGlobalMappingCreate(comm,1,maxlocal,ltog,PETSC_OWN_POINTER,mapping);
200:   PetscFree(globals);
201:   return(0);
202: }

206: /*@
207:     ISLocalToGlobalMappingGetBlockSize - Gets the blocksize of the mapping
208:     ordering and a global parallel ordering.

210:     Not Collective

212:     Input Parameters:
213: .   mapping - mapping data structure

215:     Output Parameter:
216: .   bs - the blocksize

218:     Level: advanced

220:     Concepts: mapping^local to global

222: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS()
223: @*/
224: PetscErrorCode  ISLocalToGlobalMappingGetBlockSize(ISLocalToGlobalMapping mapping,PetscInt *bs)
225: {
228:   *bs = mapping->bs;
229:   return(0);
230: }

234: /*@
235:     ISLocalToGlobalMappingCreate - Creates a mapping between a local (0 to n)
236:     ordering and a global parallel ordering.

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

240:     Input Parameters:
241: +   comm - MPI communicator
242: .   bs - the block size
243: .   n - the number of local elements divided by the block size, or equivalently the number of block indices
244: .   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
245: -   mode - see PetscCopyMode

247:     Output Parameter:
248: .   mapping - new mapping data structure

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

253:     Concepts: mapping^local to global

255: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS()
256: @*/
257: PetscErrorCode  ISLocalToGlobalMappingCreate(MPI_Comm cm,PetscInt bs,PetscInt n,const PetscInt indices[],PetscCopyMode mode,ISLocalToGlobalMapping *mapping)
258: {
260:   PetscInt       *in;


266:   *mapping = NULL;
267:   ISInitializePackage();

269:   PetscHeaderCreate(*mapping,IS_LTOGM_CLASSID,"ISLocalToGlobalMapping","Local to global mapping","IS",
270:                            cm,ISLocalToGlobalMappingDestroy,ISLocalToGlobalMappingView);
271:   (*mapping)->n             = n;
272:   (*mapping)->bs            = bs;
273:   (*mapping)->info_cached   = PETSC_FALSE;
274:   (*mapping)->info_free     = PETSC_FALSE;
275:   (*mapping)->info_procs    = NULL;
276:   (*mapping)->info_numprocs = NULL;
277:   (*mapping)->info_indices  = NULL;
278:   /*
279:     Do not create the global to local mapping. This is only created if
280:     ISGlobalToLocalMapping() is called
281:   */
282:   (*mapping)->globals = 0;
283:   if (mode == PETSC_COPY_VALUES) {
284:     PetscMalloc1(n,&in);
285:     PetscMemcpy(in,indices,n*sizeof(PetscInt));
286:     (*mapping)->indices = in;
287:     PetscLogObjectMemory((PetscObject)*mapping,n*sizeof(PetscInt));
288:   } else if (mode == PETSC_OWN_POINTER) {
289:     (*mapping)->indices = (PetscInt*)indices;
290:     PetscLogObjectMemory((PetscObject)*mapping,n*sizeof(PetscInt));
291:   }
292:   else SETERRQ(cm,PETSC_ERR_SUP,"Cannot currently use PETSC_USE_POINTER");
293:   PetscStrallocpy("basic",&((PetscObject)*mapping)->type_name);
294:   return(0);
295: }

299: /*@
300:    ISLocalToGlobalMappingDestroy - Destroys a mapping between a local (0 to n)
301:    ordering and a global parallel ordering.

303:    Note Collective

305:    Input Parameters:
306: .  mapping - mapping data structure

308:    Level: advanced

310: .seealso: ISLocalToGlobalMappingCreate()
311: @*/
312: PetscErrorCode  ISLocalToGlobalMappingDestroy(ISLocalToGlobalMapping *mapping)
313: {

317:   if (!*mapping) return(0);
319:   if (--((PetscObject)(*mapping))->refct > 0) {*mapping = 0;return(0);}
320:   PetscFree((*mapping)->indices);
321:   PetscFree((*mapping)->globals);
322:   PetscFree((*mapping)->info_procs);
323:   PetscFree((*mapping)->info_numprocs);
324:   if ((*mapping)->info_indices) {
325:     PetscInt i;

327:     PetscFree(((*mapping)->info_indices)[0]);
328:     for (i=1; i<(*mapping)->info_nproc; i++) {
329:       PetscFree(((*mapping)->info_indices)[i]);
330:     }
331:     PetscFree((*mapping)->info_indices);
332:   }
333:   PetscHeaderDestroy(mapping);
334:   *mapping = 0;
335:   return(0);
336: }

340: /*@
341:     ISLocalToGlobalMappingApplyIS - Creates from an IS in the local numbering
342:     a new index set using the global numbering defined in an ISLocalToGlobalMapping
343:     context.

345:     Not collective

347:     Input Parameters:
348: +   mapping - mapping between local and global numbering
349: -   is - index set in local numbering

351:     Output Parameters:
352: .   newis - index set in global numbering

354:     Level: advanced

356:     Concepts: mapping^local to global

358: .seealso: ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingCreate(),
359:           ISLocalToGlobalMappingDestroy(), ISGlobalToLocalMappingApply()
360: @*/
361: PetscErrorCode  ISLocalToGlobalMappingApplyIS(ISLocalToGlobalMapping mapping,IS is,IS *newis)
362: {
364:   PetscInt       n,*idxout;
365:   const PetscInt *idxin;


372:   ISGetLocalSize(is,&n);
373:   ISGetIndices(is,&idxin);
374:   PetscMalloc1(n,&idxout);
375:   ISLocalToGlobalMappingApply(mapping,n,idxin,idxout);
376:   ISRestoreIndices(is,&idxin);
377:   ISCreateGeneral(PetscObjectComm((PetscObject)is),n,idxout,PETSC_OWN_POINTER,newis);
378:   return(0);
379: }

383: /*@
384:    ISLocalToGlobalMappingApply - Takes a list of integers in a local numbering
385:    and converts them to the global numbering.

387:    Not collective

389:    Input Parameters:
390: +  mapping - the local to global mapping context
391: .  N - number of integers
392: -  in - input indices in local numbering

394:    Output Parameter:
395: .  out - indices in global numbering

397:    Notes:
398:    The in and out array parameters may be identical.

400:    Level: advanced

402: .seealso: ISLocalToGlobalMappingApplyBlock(), ISLocalToGlobalMappingCreate(),ISLocalToGlobalMappingDestroy(),
403:           ISLocalToGlobalMappingApplyIS(),AOCreateBasic(),AOApplicationToPetsc(),
404:           AOPetscToApplication(), ISGlobalToLocalMappingApply()

406:     Concepts: mapping^local to global
407: @*/
408: PetscErrorCode ISLocalToGlobalMappingApply(ISLocalToGlobalMapping mapping,PetscInt N,const PetscInt in[],PetscInt out[])
409: {
410:   PetscInt i,bs,Nmax;

414:   bs   = mapping->bs;
415:   Nmax = bs*mapping->n;
416:   if (bs == 1) {
417:     const PetscInt *idx = mapping->indices;
418:     for (i=0; i<N; i++) {
419:       if (in[i] < 0) {
420:         out[i] = in[i];
421:         continue;
422:       }
423:       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);
424:       out[i] = idx[in[i]];
425:     }
426:   } else {
427:     const PetscInt *idx = mapping->indices;
428:     for (i=0; i<N; i++) {
429:       if (in[i] < 0) {
430:         out[i] = in[i];
431:         continue;
432:       }
433:       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);
434:       out[i] = idx[in[i]/bs]*bs + (in[i] % bs);
435:     }
436:   }
437:   return(0);
438: }

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

445:    Not collective

447:    Input Parameters:
448: +  mapping - the local to global mapping context
449: .  N - number of integers
450: -  in - input indices in local block numbering

452:    Output Parameter:
453: .  out - indices in global block numbering

455:    Notes:
456:    The in and out array parameters may be identical.

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

462:    Level: advanced

464: .seealso: ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingCreate(),ISLocalToGlobalMappingDestroy(),
465:           ISLocalToGlobalMappingApplyIS(),AOCreateBasic(),AOApplicationToPetsc(),
466:           AOPetscToApplication(), ISGlobalToLocalMappingApply()

468:     Concepts: mapping^local to global
469: @*/
470: PetscErrorCode ISLocalToGlobalMappingApplyBlock(ISLocalToGlobalMapping mapping,PetscInt N,const PetscInt in[],PetscInt out[])
471: {

475:   {
476:     PetscInt i,Nmax = mapping->n;
477:     const PetscInt *idx = mapping->indices;

479:     for (i=0; i<N; i++) {
480:       if (in[i] < 0) {
481:         out[i] = in[i];
482:         continue;
483:       }
484:       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);
485:       out[i] = idx[in[i]];
486:     }
487:   }
488:   return(0);
489: }

491: /* -----------------------------------------------------------------------------------------*/

495: /*
496:     Creates the global fields in the ISLocalToGlobalMapping structure
497: */
498: static PetscErrorCode ISGlobalToLocalMappingSetUp_Private(ISLocalToGlobalMapping mapping)
499: {
501:   PetscInt       i,*idx = mapping->indices,n = mapping->n,end,start,*globals;

504:   end   = 0;
505:   start = PETSC_MAX_INT;

507:   for (i=0; i<n; i++) {
508:     if (idx[i] < 0) continue;
509:     if (idx[i] < start) start = idx[i];
510:     if (idx[i] > end)   end   = idx[i];
511:   }
512:   if (start > end) {start = 0; end = -1;}
513:   mapping->globalstart = start;
514:   mapping->globalend   = end;

516:   PetscMalloc1(end-start+2,&globals);
517:   mapping->globals = globals;
518:   for (i=0; i<end-start+1; i++) globals[i] = -1;
519:   for (i=0; i<n; i++) {
520:     if (idx[i] < 0) continue;
521:     globals[idx[i] - start] = i;
522:   }

524:   PetscLogObjectMemory((PetscObject)mapping,(end-start+1)*sizeof(PetscInt));
525:   return(0);
526: }

530: /*@
531:     ISGlobalToLocalMappingApply - Provides the local numbering for a list of integers
532:     specified with a global numbering.

534:     Not collective

536:     Input Parameters:
537: +   mapping - mapping between local and global numbering
538: .   type - IS_GTOLM_MASK - replaces global indices with no local value with -1
539:            IS_GTOLM_DROP - drops the indices with no local value from the output list
540: .   n - number of global indices to map
541: -   idx - global indices to map

543:     Output Parameters:
544: +   nout - number of indices in output array (if type == IS_GTOLM_MASK then nout = n)
545: -   idxout - local index of each global index, one must pass in an array long enough
546:              to hold all the indices. You can call ISGlobalToLocalMappingApply() with
547:              idxout == NULL to determine the required length (returned in nout)
548:              and then allocate the required space and call ISGlobalToLocalMappingApply()
549:              a second time to set the values.

551:     Notes:
552:     Either nout or idxout may be NULL. idx and idxout may be identical.

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

557:     Level: advanced

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

562:     Concepts: mapping^global to local

564: .seealso: ISLocalToGlobalMappingApply(), ISGlobalToLocalMappingApplyBlock(), ISLocalToGlobalMappingCreate(),
565:           ISLocalToGlobalMappingDestroy()
566: @*/
567: PetscErrorCode  ISGlobalToLocalMappingApply(ISLocalToGlobalMapping mapping,ISGlobalToLocalMappingType type,
568:                                             PetscInt n,const PetscInt idx[],PetscInt *nout,PetscInt idxout[])
569: {
570:   PetscInt       i,*globals,nf = 0,tmp,start,end,bs;

575:   if (!mapping->globals) {
576:     ISGlobalToLocalMappingSetUp_Private(mapping);
577:   }
578:   globals = mapping->globals;
579:   start   = mapping->globalstart;
580:   end     = mapping->globalend;
581:   bs      = mapping->bs;

583:   if (type == IS_GTOLM_MASK) {
584:     if (idxout) {
585:       for (i=0; i<n; i++) {
586:         if (idx[i] < 0)                   idxout[i] = idx[i];
587:         else if (idx[i] < bs*start)       idxout[i] = -1;
588:         else if (idx[i] > bs*(end+1)-1)   idxout[i] = -1;
589:         else                              idxout[i] = bs*globals[idx[i]/bs - start] + (idx[i] % bs);
590:       }
591:     }
592:     if (nout) *nout = n;
593:   } else {
594:     if (idxout) {
595:       for (i=0; i<n; i++) {
596:         if (idx[i] < 0) continue;
597:         if (idx[i] < bs*start) continue;
598:         if (idx[i] > bs*(end+1)-1) continue;
599:         tmp = bs*globals[idx[i]/bs - start] + (idx[i] % bs);
600:         if (tmp < 0) continue;
601:         idxout[nf++] = tmp;
602:       }
603:     } else {
604:       for (i=0; i<n; i++) {
605:         if (idx[i] < 0) continue;
606:         if (idx[i] < bs*start) continue;
607:         if (idx[i] > bs*(end+1)-1) continue;
608:         tmp = bs*globals[idx[i]/bs - start] + (idx[i] % bs);
609:         if (tmp < 0) continue;
610:         nf++;
611:       }
612:     }
613:     if (nout) *nout = nf;
614:   }
615:   return(0);
616: }

620: /*@
621:     ISGlobalToLocalMappingApplyIS - Creates from an IS in the global numbering
622:     a new index set using the local numbering defined in an ISLocalToGlobalMapping
623:     context.

625:     Not collective

627:     Input Parameters:
628: +   mapping - mapping between local and global numbering
629: -   is - index set in global numbering

631:     Output Parameters:
632: .   newis - index set in local numbering

634:     Level: advanced

636:     Concepts: mapping^local to global

638: .seealso: ISGlobalToLocalMappingApply(), ISLocalToGlobalMappingCreate(),
639:           ISLocalToGlobalMappingDestroy()
640: @*/
641: PetscErrorCode  ISGlobalToLocalMappingApplyIS(ISLocalToGlobalMapping mapping,ISGlobalToLocalMappingType type, IS is,IS *newis)
642: {
644:   PetscInt       n,nout,*idxout;
645:   const PetscInt *idxin;


652:   ISGetLocalSize(is,&n);
653:   ISGetIndices(is,&idxin);
654:   if (type == IS_GTOLM_MASK) {
655:     PetscMalloc1(n,&idxout);
656:   } else {
657:     ISGlobalToLocalMappingApply(mapping,type,n,idxin,&nout,NULL);
658:     PetscMalloc1(nout,&idxout);
659:   }
660:   ISGlobalToLocalMappingApply(mapping,type,n,idxin,&nout,idxout);
661:   ISRestoreIndices(is,&idxin);
662:   ISCreateGeneral(PetscObjectComm((PetscObject)is),nout,idxout,PETSC_OWN_POINTER,newis);
663:   return(0);
664: }

668: /*@
669:     ISGlobalToLocalMappingApplyBlock - Provides the local block numbering for a list of integers
670:     specified with a block global numbering.

672:     Not collective

674:     Input Parameters:
675: +   mapping - mapping between local and global numbering
676: .   type - IS_GTOLM_MASK - replaces global indices with no local value with -1
677:            IS_GTOLM_DROP - drops the indices with no local value from the output list
678: .   n - number of global indices to map
679: -   idx - global indices to map

681:     Output Parameters:
682: +   nout - number of indices in output array (if type == IS_GTOLM_MASK then nout = n)
683: -   idxout - local index of each global index, one must pass in an array long enough
684:              to hold all the indices. You can call ISGlobalToLocalMappingApplyBlock() with
685:              idxout == NULL to determine the required length (returned in nout)
686:              and then allocate the required space and call ISGlobalToLocalMappingApplyBlock()
687:              a second time to set the values.

689:     Notes:
690:     Either nout or idxout may be NULL. idx and idxout may be identical.

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

695:     Level: advanced

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

700:     Concepts: mapping^global to local

702: .seealso: ISLocalToGlobalMappingApply(), ISGlobalToLocalMappingApply(), ISLocalToGlobalMappingCreate(),
703:           ISLocalToGlobalMappingDestroy()
704: @*/
705: PetscErrorCode  ISGlobalToLocalMappingApplyBlock(ISLocalToGlobalMapping mapping,ISGlobalToLocalMappingType type,
706:                                   PetscInt n,const PetscInt idx[],PetscInt *nout,PetscInt idxout[])
707: {
708:   PetscInt       i,*globals,nf = 0,tmp,start,end;

713:   if (!mapping->globals) {
714:     ISGlobalToLocalMappingSetUp_Private(mapping);
715:   }
716:   globals = mapping->globals;
717:   start   = mapping->globalstart;
718:   end     = mapping->globalend;

720:   if (type == IS_GTOLM_MASK) {
721:     if (idxout) {
722:       for (i=0; i<n; i++) {
723:         if (idx[i] < 0) idxout[i] = idx[i];
724:         else if (idx[i] < start) idxout[i] = -1;
725:         else if (idx[i] > end)   idxout[i] = -1;
726:         else                     idxout[i] = globals[idx[i] - start];
727:       }
728:     }
729:     if (nout) *nout = n;
730:   } else {
731:     if (idxout) {
732:       for (i=0; i<n; i++) {
733:         if (idx[i] < 0) continue;
734:         if (idx[i] < start) continue;
735:         if (idx[i] > end) continue;
736:         tmp = globals[idx[i] - start];
737:         if (tmp < 0) continue;
738:         idxout[nf++] = tmp;
739:       }
740:     } else {
741:       for (i=0; i<n; i++) {
742:         if (idx[i] < 0) continue;
743:         if (idx[i] < start) continue;
744:         if (idx[i] > end) continue;
745:         tmp = globals[idx[i] - start];
746:         if (tmp < 0) continue;
747:         nf++;
748:       }
749:     }
750:     if (nout) *nout = nf;
751:   }
752:   return(0);
753: }

757: /*@C
758:     ISLocalToGlobalMappingGetBlockInfo - Gets the neighbor information for each processor and
759:      each index shared by more than one processor

761:     Collective on ISLocalToGlobalMapping

763:     Input Parameters:
764: .   mapping - the mapping from local to global indexing

766:     Output Parameter:
767: +   nproc - number of processors that are connected to this one
768: .   proc - neighboring processors
769: .   numproc - number of indices for each subdomain (processor)
770: -   indices - indices of nodes (in local numbering) shared with neighbors (sorted by global numbering)

772:     Level: advanced

774:     Concepts: mapping^local to global

776:     Fortran Usage:
777: $        ISLocalToGlobalMpngGetInfoSize(ISLocalToGlobalMapping,PetscInt nproc,PetscInt numprocmax,ierr) followed by
778: $        ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping,PetscInt nproc, PetscInt procs[nproc],PetscInt numprocs[nproc],
779:           PetscInt indices[nproc][numprocmax],ierr)
780:         There is no ISLocalToGlobalMappingRestoreInfo() in Fortran. You must make sure that procs[], numprocs[] and
781:         indices[][] are large enough arrays, either by allocating them dynamically or defining static ones large enough.


784: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
785:           ISLocalToGlobalMappingRestoreInfo()
786: @*/
787: PetscErrorCode  ISLocalToGlobalMappingGetBlockInfo(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
788: {

793:   if (mapping->info_cached) {
794:     *nproc = mapping->info_nproc;
795:     *procs = mapping->info_procs;
796:     *numprocs = mapping->info_numprocs;
797:     *indices = mapping->info_indices;
798:   } else {
799:     ISLocalToGlobalMappingGetBlockInfo_Private(mapping,nproc,procs,numprocs,indices);
800:   }
801:   return(0);
802: }

806: static PetscErrorCode  ISLocalToGlobalMappingGetBlockInfo_Private(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
807: {
809:   PetscMPIInt    size,rank,tag1,tag2,tag3,*len,*source,imdex;
810:   PetscInt       i,n = mapping->n,Ng,ng,max = 0,*lindices = mapping->indices;
811:   PetscInt       *nprocs,*owner,nsends,*sends,j,*starts,nmax,nrecvs,*recvs,proc;
812:   PetscInt       cnt,scale,*ownedsenders,*nownedsenders,rstart,nowned;
813:   PetscInt       node,nownedm,nt,*sends2,nsends2,*starts2,*lens2,*dest,nrecvs2,*starts3,*recvs2,k,*bprocs,*tmp;
814:   PetscInt       first_procs,first_numprocs,*first_indices;
815:   MPI_Request    *recv_waits,*send_waits;
816:   MPI_Status     recv_status,*send_status,*recv_statuses;
817:   MPI_Comm       comm;
818:   PetscBool      debug = PETSC_FALSE;

821:   PetscObjectGetComm((PetscObject)mapping,&comm);
822:   MPI_Comm_size(comm,&size);
823:   MPI_Comm_rank(comm,&rank);
824:   if (size == 1) {
825:     *nproc         = 0;
826:     *procs         = NULL;
827:     PetscMalloc(sizeof(PetscInt),numprocs);
828:     (*numprocs)[0] = 0;
829:     PetscMalloc(sizeof(PetscInt*),indices);
830:     (*indices)[0]  = NULL;
831:     /* save info for reuse */
832:     mapping->info_nproc = *nproc;
833:     mapping->info_procs = *procs;
834:     mapping->info_numprocs = *numprocs;
835:     mapping->info_indices = *indices;
836:     mapping->info_cached = PETSC_TRUE;
837:     return(0);
838:   }

840:   PetscOptionsGetBool(NULL,"-islocaltoglobalmappinggetinfo_debug",&debug,NULL);

842:   /*
843:     Notes on ISLocalToGlobalMappingGetBlockInfo

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

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

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

853:     nontrivial locally owned node - node that is not in the interior (i.e. has more than one
854:            local subdomain
855:   */
856:   PetscObjectGetNewTag((PetscObject)mapping,&tag1);
857:   PetscObjectGetNewTag((PetscObject)mapping,&tag2);
858:   PetscObjectGetNewTag((PetscObject)mapping,&tag3);

860:   for (i=0; i<n; i++) {
861:     if (lindices[i] > max) max = lindices[i];
862:   }
863:   MPI_Allreduce(&max,&Ng,1,MPIU_INT,MPI_MAX,comm);
864:   Ng++;
865:   MPI_Comm_size(comm,&size);
866:   MPI_Comm_rank(comm,&rank);
867:   scale  = Ng/size + 1;
868:   ng     = scale; if (rank == size-1) ng = Ng - scale*(size-1); ng = PetscMax(1,ng);
869:   rstart = scale*rank;

871:   /* determine ownership ranges of global indices */
872:   PetscMalloc1(2*size,&nprocs);
873:   PetscMemzero(nprocs,2*size*sizeof(PetscInt));

875:   /* determine owners of each local node  */
876:   PetscMalloc1(n,&owner);
877:   for (i=0; i<n; i++) {
878:     proc             = lindices[i]/scale; /* processor that globally owns this index */
879:     nprocs[2*proc+1] = 1;                 /* processor globally owns at least one of ours */
880:     owner[i]         = proc;
881:     nprocs[2*proc]++;                     /* count of how many that processor globally owns of ours */
882:   }
883:   nsends = 0; for (i=0; i<size; i++) nsends += nprocs[2*i+1];
884:   PetscInfo1(mapping,"Number of global owners for my local data %D\n",nsends);

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

890:   /* post receives for owned rows */
891:   PetscMalloc1((2*nrecvs+1)*(nmax+1),&recvs);
892:   PetscMalloc1(nrecvs+1,&recv_waits);
893:   for (i=0; i<nrecvs; i++) {
894:     MPI_Irecv(recvs+2*nmax*i,2*nmax,MPIU_INT,MPI_ANY_SOURCE,tag1,comm,recv_waits+i);
895:   }

897:   /* pack messages containing lists of local nodes to owners */
898:   PetscMalloc1(2*n+1,&sends);
899:   PetscMalloc1(size+1,&starts);
900:   starts[0] = 0;
901:   for (i=1; i<size; i++) starts[i] = starts[i-1] + 2*nprocs[2*i-2];
902:   for (i=0; i<n; i++) {
903:     sends[starts[owner[i]]++] = lindices[i];
904:     sends[starts[owner[i]]++] = i;
905:   }
906:   PetscFree(owner);
907:   starts[0] = 0;
908:   for (i=1; i<size; i++) starts[i] = starts[i-1] + 2*nprocs[2*i-2];

910:   /* send the messages */
911:   PetscMalloc1(nsends+1,&send_waits);
912:   PetscMalloc1(nsends+1,&dest);
913:   cnt = 0;
914:   for (i=0; i<size; i++) {
915:     if (nprocs[2*i]) {
916:       MPI_Isend(sends+starts[i],2*nprocs[2*i],MPIU_INT,i,tag1,comm,send_waits+cnt);
917:       dest[cnt] = i;
918:       cnt++;
919:     }
920:   }
921:   PetscFree(starts);

923:   /* wait on receives */
924:   PetscMalloc1(nrecvs+1,&source);
925:   PetscMalloc1(nrecvs+1,&len);
926:   cnt  = nrecvs;
927:   PetscMalloc1(ng+1,&nownedsenders);
928:   PetscMemzero(nownedsenders,ng*sizeof(PetscInt));
929:   while (cnt) {
930:     MPI_Waitany(nrecvs,recv_waits,&imdex,&recv_status);
931:     /* unpack receives into our local space */
932:     MPI_Get_count(&recv_status,MPIU_INT,&len[imdex]);
933:     source[imdex] = recv_status.MPI_SOURCE;
934:     len[imdex]    = len[imdex]/2;
935:     /* count how many local owners for each of my global owned indices */
936:     for (i=0; i<len[imdex]; i++) nownedsenders[recvs[2*imdex*nmax+2*i]-rstart]++;
937:     cnt--;
938:   }
939:   PetscFree(recv_waits);

941:   /* count how many globally owned indices are on an edge multiplied by how many processors own them. */
942:   nowned  = 0;
943:   nownedm = 0;
944:   for (i=0; i<ng; i++) {
945:     if (nownedsenders[i] > 1) {nownedm += nownedsenders[i]; nowned++;}
946:   }

948:   /* create single array to contain rank of all local owners of each globally owned index */
949:   PetscMalloc1(nownedm+1,&ownedsenders);
950:   PetscMalloc1(ng+1,&starts);
951:   starts[0] = 0;
952:   for (i=1; i<ng; i++) {
953:     if (nownedsenders[i-1] > 1) starts[i] = starts[i-1] + nownedsenders[i-1];
954:     else starts[i] = starts[i-1];
955:   }

957:   /* for each nontrival globally owned node list all arriving processors */
958:   for (i=0; i<nrecvs; i++) {
959:     for (j=0; j<len[i]; j++) {
960:       node = recvs[2*i*nmax+2*j]-rstart;
961:       if (nownedsenders[node] > 1) ownedsenders[starts[node]++] = source[i];
962:     }
963:   }

965:   if (debug) { /* -----------------------------------  */
966:     starts[0] = 0;
967:     for (i=1; i<ng; i++) {
968:       if (nownedsenders[i-1] > 1) starts[i] = starts[i-1] + nownedsenders[i-1];
969:       else starts[i] = starts[i-1];
970:     }
971:     for (i=0; i<ng; i++) {
972:       if (nownedsenders[i] > 1) {
973:         PetscSynchronizedPrintf(comm,"[%d] global node %D local owner processors: ",rank,i+rstart);
974:         for (j=0; j<nownedsenders[i]; j++) {
975:           PetscSynchronizedPrintf(comm,"%D ",ownedsenders[starts[i]+j]);
976:         }
977:         PetscSynchronizedPrintf(comm,"\n");
978:       }
979:     }
980:     PetscSynchronizedFlush(comm,PETSC_STDOUT);
981:   } /* -----------------------------------  */

983:   /* wait on original sends */
984:   if (nsends) {
985:     PetscMalloc1(nsends,&send_status);
986:     MPI_Waitall(nsends,send_waits,send_status);
987:     PetscFree(send_status);
988:   }
989:   PetscFree(send_waits);
990:   PetscFree(sends);
991:   PetscFree(nprocs);

993:   /* pack messages to send back to local owners */
994:   starts[0] = 0;
995:   for (i=1; i<ng; i++) {
996:     if (nownedsenders[i-1] > 1) starts[i] = starts[i-1] + nownedsenders[i-1];
997:     else starts[i] = starts[i-1];
998:   }
999:   nsends2 = nrecvs;
1000:   PetscMalloc1(nsends2+1,&nprocs); /* length of each message */
1001:   for (i=0; i<nrecvs; i++) {
1002:     nprocs[i] = 1;
1003:     for (j=0; j<len[i]; j++) {
1004:       node = recvs[2*i*nmax+2*j]-rstart;
1005:       if (nownedsenders[node] > 1) nprocs[i] += 2 + nownedsenders[node];
1006:     }
1007:   }
1008:   nt = 0;
1009:   for (i=0; i<nsends2; i++) nt += nprocs[i];

1011:   PetscMalloc1(nt+1,&sends2);
1012:   PetscMalloc1(nsends2+1,&starts2);

1014:   starts2[0] = 0;
1015:   for (i=1; i<nsends2; i++) starts2[i] = starts2[i-1] + nprocs[i-1];
1016:   /*
1017:      Each message is 1 + nprocs[i] long, and consists of
1018:        (0) the number of nodes being sent back
1019:        (1) the local node number,
1020:        (2) the number of processors sharing it,
1021:        (3) the processors sharing it
1022:   */
1023:   for (i=0; i<nsends2; i++) {
1024:     cnt = 1;
1025:     sends2[starts2[i]] = 0;
1026:     for (j=0; j<len[i]; j++) {
1027:       node = recvs[2*i*nmax+2*j]-rstart;
1028:       if (nownedsenders[node] > 1) {
1029:         sends2[starts2[i]]++;
1030:         sends2[starts2[i]+cnt++] = recvs[2*i*nmax+2*j+1];
1031:         sends2[starts2[i]+cnt++] = nownedsenders[node];
1032:         PetscMemcpy(&sends2[starts2[i]+cnt],&ownedsenders[starts[node]],nownedsenders[node]*sizeof(PetscInt));
1033:         cnt += nownedsenders[node];
1034:       }
1035:     }
1036:   }

1038:   /* receive the message lengths */
1039:   nrecvs2 = nsends;
1040:   PetscMalloc1(nrecvs2+1,&lens2);
1041:   PetscMalloc1(nrecvs2+1,&starts3);
1042:   PetscMalloc1(nrecvs2+1,&recv_waits);
1043:   for (i=0; i<nrecvs2; i++) {
1044:     MPI_Irecv(&lens2[i],1,MPIU_INT,dest[i],tag2,comm,recv_waits+i);
1045:   }

1047:   /* send the message lengths */
1048:   for (i=0; i<nsends2; i++) {
1049:     MPI_Send(&nprocs[i],1,MPIU_INT,source[i],tag2,comm);
1050:   }

1052:   /* wait on receives of lens */
1053:   if (nrecvs2) {
1054:     PetscMalloc1(nrecvs2,&recv_statuses);
1055:     MPI_Waitall(nrecvs2,recv_waits,recv_statuses);
1056:     PetscFree(recv_statuses);
1057:   }
1058:   PetscFree(recv_waits);

1060:   starts3[0] = 0;
1061:   nt         = 0;
1062:   for (i=0; i<nrecvs2-1; i++) {
1063:     starts3[i+1] = starts3[i] + lens2[i];
1064:     nt          += lens2[i];
1065:   }
1066:   if (nrecvs2) nt += lens2[nrecvs2-1];

1068:   PetscMalloc1(nt+1,&recvs2);
1069:   PetscMalloc1(nrecvs2+1,&recv_waits);
1070:   for (i=0; i<nrecvs2; i++) {
1071:     MPI_Irecv(recvs2+starts3[i],lens2[i],MPIU_INT,dest[i],tag3,comm,recv_waits+i);
1072:   }

1074:   /* send the messages */
1075:   PetscMalloc1(nsends2+1,&send_waits);
1076:   for (i=0; i<nsends2; i++) {
1077:     MPI_Isend(sends2+starts2[i],nprocs[i],MPIU_INT,source[i],tag3,comm,send_waits+i);
1078:   }

1080:   /* wait on receives */
1081:   if (nrecvs2) {
1082:     PetscMalloc1(nrecvs2,&recv_statuses);
1083:     MPI_Waitall(nrecvs2,recv_waits,recv_statuses);
1084:     PetscFree(recv_statuses);
1085:   }
1086:   PetscFree(recv_waits);
1087:   PetscFree(nprocs);

1089:   if (debug) { /* -----------------------------------  */
1090:     cnt = 0;
1091:     for (i=0; i<nrecvs2; i++) {
1092:       nt = recvs2[cnt++];
1093:       for (j=0; j<nt; j++) {
1094:         PetscSynchronizedPrintf(comm,"[%d] local node %D number of subdomains %D: ",rank,recvs2[cnt],recvs2[cnt+1]);
1095:         for (k=0; k<recvs2[cnt+1]; k++) {
1096:           PetscSynchronizedPrintf(comm,"%D ",recvs2[cnt+2+k]);
1097:         }
1098:         cnt += 2 + recvs2[cnt+1];
1099:         PetscSynchronizedPrintf(comm,"\n");
1100:       }
1101:     }
1102:     PetscSynchronizedFlush(comm,PETSC_STDOUT);
1103:   } /* -----------------------------------  */

1105:   /* count number subdomains for each local node */
1106:   PetscMalloc1(size,&nprocs);
1107:   PetscMemzero(nprocs,size*sizeof(PetscInt));
1108:   cnt  = 0;
1109:   for (i=0; i<nrecvs2; i++) {
1110:     nt = recvs2[cnt++];
1111:     for (j=0; j<nt; j++) {
1112:       for (k=0; k<recvs2[cnt+1]; k++) nprocs[recvs2[cnt+2+k]]++;
1113:       cnt += 2 + recvs2[cnt+1];
1114:     }
1115:   }
1116:   nt = 0; for (i=0; i<size; i++) nt += (nprocs[i] > 0);
1117:   *nproc    = nt;
1118:   PetscMalloc1(nt+1,procs);
1119:   PetscMalloc1(nt+1,numprocs);
1120:   PetscMalloc1(nt+1,indices);
1121:   for (i=0;i<nt+1;i++) (*indices)[i]=NULL;
1122:   PetscMalloc1(size,&bprocs);
1123:   cnt  = 0;
1124:   for (i=0; i<size; i++) {
1125:     if (nprocs[i] > 0) {
1126:       bprocs[i]        = cnt;
1127:       (*procs)[cnt]    = i;
1128:       (*numprocs)[cnt] = nprocs[i];
1129:       PetscMalloc1(nprocs[i],&(*indices)[cnt]);
1130:       cnt++;
1131:     }
1132:   }

1134:   /* make the list of subdomains for each nontrivial local node */
1135:   PetscMemzero(*numprocs,nt*sizeof(PetscInt));
1136:   cnt  = 0;
1137:   for (i=0; i<nrecvs2; i++) {
1138:     nt = recvs2[cnt++];
1139:     for (j=0; j<nt; j++) {
1140:       for (k=0; k<recvs2[cnt+1]; k++) (*indices)[bprocs[recvs2[cnt+2+k]]][(*numprocs)[bprocs[recvs2[cnt+2+k]]]++] = recvs2[cnt];
1141:       cnt += 2 + recvs2[cnt+1];
1142:     }
1143:   }
1144:   PetscFree(bprocs);
1145:   PetscFree(recvs2);

1147:   /* sort the node indexing by their global numbers */
1148:   nt = *nproc;
1149:   for (i=0; i<nt; i++) {
1150:     PetscMalloc1((*numprocs)[i],&tmp);
1151:     for (j=0; j<(*numprocs)[i]; j++) tmp[j] = lindices[(*indices)[i][j]];
1152:     PetscSortIntWithArray((*numprocs)[i],tmp,(*indices)[i]);
1153:     PetscFree(tmp);
1154:   }

1156:   if (debug) { /* -----------------------------------  */
1157:     nt = *nproc;
1158:     for (i=0; i<nt; i++) {
1159:       PetscSynchronizedPrintf(comm,"[%d] subdomain %D number of indices %D: ",rank,(*procs)[i],(*numprocs)[i]);
1160:       for (j=0; j<(*numprocs)[i]; j++) {
1161:         PetscSynchronizedPrintf(comm,"%D ",(*indices)[i][j]);
1162:       }
1163:       PetscSynchronizedPrintf(comm,"\n");
1164:     }
1165:     PetscSynchronizedFlush(comm,PETSC_STDOUT);
1166:   } /* -----------------------------------  */

1168:   /* wait on sends */
1169:   if (nsends2) {
1170:     PetscMalloc1(nsends2,&send_status);
1171:     MPI_Waitall(nsends2,send_waits,send_status);
1172:     PetscFree(send_status);
1173:   }

1175:   PetscFree(starts3);
1176:   PetscFree(dest);
1177:   PetscFree(send_waits);

1179:   PetscFree(nownedsenders);
1180:   PetscFree(ownedsenders);
1181:   PetscFree(starts);
1182:   PetscFree(starts2);
1183:   PetscFree(lens2);

1185:   PetscFree(source);
1186:   PetscFree(len);
1187:   PetscFree(recvs);
1188:   PetscFree(nprocs);
1189:   PetscFree(sends2);

1191:   /* put the information about myself as the first entry in the list */
1192:   first_procs    = (*procs)[0];
1193:   first_numprocs = (*numprocs)[0];
1194:   first_indices  = (*indices)[0];
1195:   for (i=0; i<*nproc; i++) {
1196:     if ((*procs)[i] == rank) {
1197:       (*procs)[0]    = (*procs)[i];
1198:       (*numprocs)[0] = (*numprocs)[i];
1199:       (*indices)[0]  = (*indices)[i];
1200:       (*procs)[i]    = first_procs;
1201:       (*numprocs)[i] = first_numprocs;
1202:       (*indices)[i]  = first_indices;
1203:       break;
1204:     }
1205:   }

1207:   /* save info for reuse */
1208:   mapping->info_nproc = *nproc;
1209:   mapping->info_procs = *procs;
1210:   mapping->info_numprocs = *numprocs;
1211:   mapping->info_indices = *indices;
1212:   mapping->info_cached = PETSC_TRUE;
1213:   return(0);
1214: }

1218: /*@C
1219:     ISLocalToGlobalMappingRestoreBlockInfo - Frees the memory allocated by ISLocalToGlobalMappingGetBlockInfo()

1221:     Collective on ISLocalToGlobalMapping

1223:     Input Parameters:
1224: .   mapping - the mapping from local to global indexing

1226:     Output Parameter:
1227: +   nproc - number of processors that are connected to this one
1228: .   proc - neighboring processors
1229: .   numproc - number of indices for each processor
1230: -   indices - indices of local nodes shared with neighbor (sorted by global numbering)

1232:     Level: advanced

1234: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
1235:           ISLocalToGlobalMappingGetInfo()
1236: @*/
1237: PetscErrorCode  ISLocalToGlobalMappingRestoreBlockInfo(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
1238: {

1243:   if (mapping->info_free) {
1244:     PetscFree(*numprocs);
1245:     if (*indices) {
1246:       PetscInt i;

1248:       PetscFree((*indices)[0]);
1249:       for (i=1; i<*nproc; i++) {
1250:         PetscFree((*indices)[i]);
1251:       }
1252:       PetscFree(*indices);
1253:     }
1254:   }
1255:   *nproc = 0;
1256:   *procs = NULL;
1257:   *numprocs = NULL;
1258:   *indices = NULL;
1259:   return(0);
1260: }

1264: /*@C
1265:     ISLocalToGlobalMappingGetInfo - Gets the neighbor information for each processor and
1266:      each index shared by more than one processor

1268:     Collective on ISLocalToGlobalMapping

1270:     Input Parameters:
1271: .   mapping - the mapping from local to global indexing

1273:     Output Parameter:
1274: +   nproc - number of processors that are connected to this one
1275: .   proc - neighboring processors
1276: .   numproc - number of indices for each subdomain (processor)
1277: -   indices - indices of nodes (in local numbering) shared with neighbors (sorted by global numbering)

1279:     Level: advanced

1281:     Concepts: mapping^local to global

1283:     Fortran Usage:
1284: $        ISLocalToGlobalMpngGetInfoSize(ISLocalToGlobalMapping,PetscInt nproc,PetscInt numprocmax,ierr) followed by
1285: $        ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping,PetscInt nproc, PetscInt procs[nproc],PetscInt numprocs[nproc],
1286:           PetscInt indices[nproc][numprocmax],ierr)
1287:         There is no ISLocalToGlobalMappingRestoreInfo() in Fortran. You must make sure that procs[], numprocs[] and
1288:         indices[][] are large enough arrays, either by allocating them dynamically or defining static ones large enough.


1291: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
1292:           ISLocalToGlobalMappingRestoreInfo()
1293: @*/
1294: PetscErrorCode  ISLocalToGlobalMappingGetInfo(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
1295: {
1297:   PetscInt       **bindices = NULL,*bnumprocs = NULL,bs = mapping->bs,i,j,k;

1301:   ISLocalToGlobalMappingGetBlockInfo(mapping,nproc,procs,&bnumprocs,&bindices);
1302:   if (bs > 1) { /* we need to expand the cached info */
1303:     PetscCalloc1(*nproc,&*indices);
1304:     PetscCalloc1(*nproc,&*numprocs);
1305:     for (i=0; i<*nproc; i++) {
1306:       PetscMalloc1(bs*bnumprocs[i],&(*indices)[i]);
1307:       for (j=0; j<bnumprocs[i]; j++) {
1308:         for (k=0; k<bs; k++) {
1309:           (*indices)[i][j*bs+k] = bs*bindices[i][j] + k;
1310:         }
1311:       }
1312:       (*numprocs)[i] = bnumprocs[i]*bs;
1313:     }
1314:     mapping->info_free = PETSC_TRUE;
1315:   } else {
1316:     *numprocs = bnumprocs;
1317:     *indices  = bindices;
1318:   }
1319:   return(0);
1320: }

1324: /*@C
1325:     ISLocalToGlobalMappingRestoreInfo - Frees the memory allocated by ISLocalToGlobalMappingGetInfo()

1327:     Collective on ISLocalToGlobalMapping

1329:     Input Parameters:
1330: .   mapping - the mapping from local to global indexing

1332:     Output Parameter:
1333: +   nproc - number of processors that are connected to this one
1334: .   proc - neighboring processors
1335: .   numproc - number of indices for each processor
1336: -   indices - indices of local nodes shared with neighbor (sorted by global numbering)

1338:     Level: advanced

1340: .seealso: ISLocalToGlobalMappingDestroy(), ISLocalToGlobalMappingCreateIS(), ISLocalToGlobalMappingCreate(),
1341:           ISLocalToGlobalMappingGetInfo()
1342: @*/
1343: PetscErrorCode  ISLocalToGlobalMappingRestoreInfo(ISLocalToGlobalMapping mapping,PetscInt *nproc,PetscInt *procs[],PetscInt *numprocs[],PetscInt **indices[])
1344: {

1348:   ISLocalToGlobalMappingRestoreBlockInfo(mapping,nproc,procs,numprocs,indices);
1349:   return(0);
1350: }

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

1357:    Not Collective

1359:    Input Arguments:
1360: . ltog - local to global mapping

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

1365:    Level: advanced

1367:    Notes: ISLocalToGlobalMappingGetSize() returns the length the this array

1369: .seealso: ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingRestoreIndices(), ISLocalToGlobalMappingGetBlockIndices(), ISLocalToGlobalMappingRestoreBlockIndices()
1370: @*/
1371: PetscErrorCode  ISLocalToGlobalMappingGetIndices(ISLocalToGlobalMapping ltog,const PetscInt **array)
1372: {
1376:   if (ltog->bs == 1) {
1377:     *array = ltog->indices;
1378:   } else {
1379:     PetscInt       *jj,k,i,j,n = ltog->n, bs = ltog->bs;
1380:     const PetscInt *ii;

1383:     PetscMalloc1(bs*n,&jj);
1384:     *array = jj;
1385:     k    = 0;
1386:     ii   = ltog->indices;
1387:     for (i=0; i<n; i++)
1388:       for (j=0; j<bs; j++)
1389:         jj[k++] = bs*ii[i] + j;
1390:   }
1391:   return(0);
1392: }

1396: /*@C
1397:    ISLocalToGlobalMappingRestoreIndices - Restore indices obtained with ISLocalToGlobalMappingRestoreIndices()

1399:    Not Collective

1401:    Input Arguments:
1402: + ltog - local to global mapping
1403: - array - array of indices

1405:    Level: advanced

1407: .seealso: ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingGetIndices()
1408: @*/
1409: PetscErrorCode  ISLocalToGlobalMappingRestoreIndices(ISLocalToGlobalMapping ltog,const PetscInt **array)
1410: {
1414:   if (ltog->bs == 1 && *array != ltog->indices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_BADPTR,"Trying to return mismatched pointer");

1416:   if (ltog->bs > 1) {
1418:     PetscFree(*(void**)array);
1419:   }
1420:   return(0);
1421: }

1425: /*@C
1426:    ISLocalToGlobalMappingGetBlockIndices - Get global indices for every local block

1428:    Not Collective

1430:    Input Arguments:
1431: . ltog - local to global mapping

1433:    Output Arguments:
1434: . array - array of indices

1436:    Level: advanced

1438: .seealso: ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingRestoreBlockIndices()
1439: @*/
1440: PetscErrorCode  ISLocalToGlobalMappingGetBlockIndices(ISLocalToGlobalMapping ltog,const PetscInt **array)
1441: {
1445:   *array = ltog->indices;
1446:   return(0);
1447: }

1451: /*@C
1452:    ISLocalToGlobalMappingRestoreBlockIndices - Restore indices obtained with ISLocalToGlobalMappingGetBlockIndices()

1454:    Not Collective

1456:    Input Arguments:
1457: + ltog - local to global mapping
1458: - array - array of indices

1460:    Level: advanced

1462: .seealso: ISLocalToGlobalMappingCreate(), ISLocalToGlobalMappingApply(), ISLocalToGlobalMappingGetIndices()
1463: @*/
1464: PetscErrorCode  ISLocalToGlobalMappingRestoreBlockIndices(ISLocalToGlobalMapping ltog,const PetscInt **array)
1465: {
1469:   if (*array != ltog->indices) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_BADPTR,"Trying to return mismatched pointer");
1470:   *array = NULL;
1471:   return(0);
1472: }

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

1479:    Not Collective

1481:    Input Arguments:
1482: + comm - communicator for the new mapping, must contain the communicator of every mapping to concatenate
1483: . n - number of mappings to concatenate
1484: - ltogs - local to global mappings

1486:    Output Arguments:
1487: . ltogcat - new mapping

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

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

1493:    Level: advanced

1495: .seealso: ISLocalToGlobalMappingCreate()
1496: @*/
1497: PetscErrorCode ISLocalToGlobalMappingConcatenate(MPI_Comm comm,PetscInt n,const ISLocalToGlobalMapping ltogs[],ISLocalToGlobalMapping *ltogcat)
1498: {
1499:   PetscInt       i,cnt,m,*idx;

1503:   if (n < 0) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have a non-negative number of mappings, given %D",n);
1507:   for (cnt=0,i=0; i<n; i++) {
1508:     ISLocalToGlobalMappingGetSize(ltogs[i],&m);
1509:     cnt += m;
1510:   }
1511:   PetscMalloc1(cnt,&idx);
1512:   for (cnt=0,i=0; i<n; i++) {
1513:     const PetscInt *subidx;
1514:     ISLocalToGlobalMappingGetSize(ltogs[i],&m);
1515:     ISLocalToGlobalMappingGetIndices(ltogs[i],&subidx);
1516:     PetscMemcpy(&idx[cnt],subidx,m*sizeof(PetscInt));
1517:     ISLocalToGlobalMappingRestoreIndices(ltogs[i],&subidx);
1518:     cnt += m;
1519:   }
1520:   ISLocalToGlobalMappingCreate(comm,1,cnt,idx,PETSC_OWN_POINTER,ltogcat);
1521:   return(0);
1522: }