Actual source code: index.c

petsc-master 2019-09-15
Report Typos and Errors

  2: /*
  3:    Defines the abstract operations on index sets, i.e. the public interface.
  4: */
  5:  #include <petsc/private/isimpl.h>
  6:  #include <petscviewer.h>
  7:  #include <petscsf.h>

  9: /* Logging support */
 10: PetscClassId IS_CLASSID;

 12: /*@
 13:    ISRenumber - Renumbers an index set (with multiplicities) in a contiguous way.

 15:    Collective on IS

 17:    Input Parmeters:
 18: +  subset - the index set
 19: -  subset_mult - the multiplcity of each entry in subset (optional, can be NULL)

 21:    Output Parameters:
 22: +  N - the maximum entry of the new IS
 23: -  subset_n - the new IS

 25:    Level: intermediate

 27: .seealso:
 28: @*/
 29: PetscErrorCode ISRenumber(IS subset, IS subset_mult, PetscInt *N, IS *subset_n)
 30: {
 31:   PetscSF        sf;
 32:   PetscLayout    map;
 33:   const PetscInt *idxs;
 34:   PetscInt       *leaf_data,*root_data,*gidxs;
 35:   PetscInt       N_n,n,i,lbounds[2],gbounds[2],Nl;
 36:   PetscInt       n_n,nlocals,start,first_index;
 37:   PetscMPIInt    commsize;
 38:   PetscBool      first_found;

 43:   if (subset_mult) {
 45:   }
 46:   if (!N && !subset_n) return(0);
 47:   ISGetLocalSize(subset,&n);
 48:   if (subset_mult) {
 49:     ISGetLocalSize(subset,&i);
 50:     if (i != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Local subset and multiplicity sizes don't match! %d != %d",n,i);
 51:   }
 52:   /* create workspace layout for computing global indices of subset */
 53:   ISGetIndices(subset,&idxs);
 54:   lbounds[0] = lbounds[1] = 0;
 55:   for (i=0;i<n;i++) {
 56:     if (idxs[i] < lbounds[0]) lbounds[0] = idxs[i];
 57:     else if (idxs[i] > lbounds[1]) lbounds[1] = idxs[i];
 58:   }
 59:   lbounds[0] = -lbounds[0];
 60:   MPIU_Allreduce(lbounds,gbounds,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)subset));
 61:   gbounds[0] = -gbounds[0];
 62:   N_n= gbounds[1] - gbounds[0] + 1;
 63:   PetscLayoutCreate(PetscObjectComm((PetscObject)subset),&map);
 64:   PetscLayoutSetBlockSize(map,1);
 65:   PetscLayoutSetSize(map,N_n);
 66:   PetscLayoutSetUp(map);
 67:   PetscLayoutGetLocalSize(map,&Nl);

 69:   /* create sf : leaf_data == multiplicity of indexes, root data == global index in layout */
 70:   PetscMalloc2(n,&leaf_data,Nl,&root_data);
 71:   if (subset_mult) {
 72:     const PetscInt* idxs_mult;

 74:     ISGetIndices(subset_mult,&idxs_mult);
 75:     PetscArraycpy(leaf_data,idxs_mult,n);
 76:     ISRestoreIndices(subset_mult,&idxs_mult);
 77:   } else {
 78:     for (i=0;i<n;i++) leaf_data[i] = 1;
 79:   }
 80:   /* local size of new subset */
 81:   n_n = 0;
 82:   for (i=0;i<n;i++) n_n += leaf_data[i];

 84:   /* global indexes in layout */
 85:   PetscMalloc1(n_n,&gidxs); /* allocating possibly extra space in gidxs which will be used later */
 86:   for (i=0;i<n;i++) gidxs[i] = idxs[i] - gbounds[0];
 87:   ISRestoreIndices(subset,&idxs);
 88:   PetscSFCreate(PetscObjectComm((PetscObject)subset),&sf);
 89:   PetscSFSetGraphLayout(sf,map,n,NULL,PETSC_COPY_VALUES,gidxs);
 90:   PetscLayoutDestroy(&map);

 92:   /* reduce from leaves to roots */
 93:   PetscArrayzero(root_data,Nl);
 94:   PetscSFReduceBegin(sf,MPIU_INT,leaf_data,root_data,MPI_MAX);
 95:   PetscSFReduceEnd(sf,MPIU_INT,leaf_data,root_data,MPI_MAX);

 97:   /* count indexes in local part of layout */
 98:   nlocals = 0;
 99:   first_index = -1;
100:   first_found = PETSC_FALSE;
101:   for (i=0;i<Nl;i++) {
102:     if (!first_found && root_data[i]) {
103:       first_found = PETSC_TRUE;
104:       first_index = i;
105:     }
106:     nlocals += root_data[i];
107:   }

109:   /* cumulative of number of indexes and size of subset without holes */
110: #if defined(PETSC_HAVE_MPI_EXSCAN)
111:   start = 0;
112:   MPI_Exscan(&nlocals,&start,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)subset));
113: #else
114:   MPI_Scan(&nlocals,&start,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)subset));
115:   start = start-nlocals;
116: #endif

118:   if (N) { /* compute total size of new subset if requested */
119:     *N   = start + nlocals;
120:     MPI_Comm_size(PetscObjectComm((PetscObject)subset),&commsize);
121:     MPI_Bcast(N,1,MPIU_INT,commsize-1,PetscObjectComm((PetscObject)subset));
122:   }

124:   if (!subset_n) {
125:     PetscFree(gidxs);
126:     PetscSFDestroy(&sf);
127:     PetscFree2(leaf_data,root_data);
128:     return(0);
129:   }

131:   /* adapt root data with cumulative */
132:   if (first_found) {
133:     PetscInt old_index;

135:     root_data[first_index] += start;
136:     old_index = first_index;
137:     for (i=first_index+1;i<Nl;i++) {
138:       if (root_data[i]) {
139:         root_data[i] += root_data[old_index];
140:         old_index = i;
141:       }
142:     }
143:   }

145:   /* from roots to leaves */
146:   PetscSFBcastBegin(sf,MPIU_INT,root_data,leaf_data);
147:   PetscSFBcastEnd(sf,MPIU_INT,root_data,leaf_data);
148:   PetscSFDestroy(&sf);

150:   /* create new IS with global indexes without holes */
151:   if (subset_mult) {
152:     const PetscInt* idxs_mult;
153:     PetscInt        cum;

155:     cum = 0;
156:     ISGetIndices(subset_mult,&idxs_mult);
157:     for (i=0;i<n;i++) {
158:       PetscInt j;
159:       for (j=0;j<idxs_mult[i];j++) gidxs[cum++] = leaf_data[i] - idxs_mult[i] + j;
160:     }
161:     ISRestoreIndices(subset_mult,&idxs_mult);
162:   } else {
163:     for (i=0;i<n;i++) {
164:       gidxs[i] = leaf_data[i]-1;
165:     }
166:   }
167:   ISCreateGeneral(PetscObjectComm((PetscObject)subset),n_n,gidxs,PETSC_OWN_POINTER,subset_n);
168:   PetscFree2(leaf_data,root_data);
169:   return(0);
170: }


173: /*@
174:    ISCreateSubIS - Create a sub index set from a global index set selecting some components.

176:    Collective on IS

178:    Input Parmeters:
179: +  is - the index set
180: -  comps - which components we will extract from is

182:    Output Parameters:
183: .  subis - the new sub index set

185:    Level: intermediate

187:    Example usage:
188:    We have an index set (is) living on 3 processes with the following values:
189:    | 4 9 0 | 2 6 7 | 10 11 1|
190:    and another index set (comps) used to indicate which components of is  we want to take,
191:    | 7 5  | 1 2 | 0 4|
192:    The output index set (subis) should look like:
193:    | 11 7 | 9 0 | 4 6|

195: .seealso: VecGetSubVector(), MatCreateSubMatrix()
196: @*/
197: PetscErrorCode ISCreateSubIS(IS is,IS comps,IS *subis)
198: {
199:   PetscSF         sf;
200:   const PetscInt  *is_indices,*comps_indices;
201:   PetscInt        *subis_indices,nroots,nleaves,*mine,i,owner,lidx;
202:   PetscSFNode     *remote;
203:   PetscErrorCode  ierr;
204:   MPI_Comm        comm;


211:   PetscObjectGetComm((PetscObject)is, &comm);
212:   ISGetLocalSize(comps,&nleaves);
213:   ISGetLocalSize(is,&nroots);
214:   PetscMalloc1(nleaves,&remote);
215:   PetscMalloc1(nleaves,&mine);
216:   ISGetIndices(comps,&comps_indices);
217:   /*
218:    * Construct a PetscSF in which "is" data serves as roots and "subis" is leaves.
219:    * Root data are sent to leaves using PetscSFBcast().
220:    * */
221:   for (i=0; i<nleaves; i++) {
222:     mine[i] = i;
223:     /* Connect a remote root with the current leaf. The value on the remote root
224:      * will be received by the current local leaf.
225:      * */
226:     owner = -1;
227:     lidx =  -1;
228:     PetscLayoutFindOwnerIndex(is->map,comps_indices[i],&owner, &lidx);
229:     remote[i].rank = owner;
230:     remote[i].index = lidx;
231:   }
232:   ISRestoreIndices(comps,&comps_indices);
233:   PetscSFCreate(comm,&sf);
234:   PetscSFSetFromOptions(sf);\
235:   PetscSFSetGraph(sf,nroots,nleaves,mine,PETSC_OWN_POINTER,remote,PETSC_OWN_POINTER);

237:   PetscMalloc1(nleaves,&subis_indices);
238:   ISGetIndices(is, &is_indices);
239:   PetscSFBcastBegin(sf,MPIU_INT,is_indices,subis_indices);
240:   PetscSFBcastEnd(sf,MPIU_INT,is_indices,subis_indices);
241:   ISRestoreIndices(is,&is_indices);
242:   PetscSFDestroy(&sf);
243:   ISCreateGeneral(comm,nleaves,subis_indices,PETSC_OWN_POINTER,subis);
244:   return(0);
245: }


248: /*@
249:    ISIdentity - Determines whether index set is the identity mapping.

251:    Collective on IS

253:    Input Parmeters:
254: .  is - the index set

256:    Output Parameters:
257: .  ident - PETSC_TRUE if an identity, else PETSC_FALSE

259:    Level: intermediate


262: .seealso: ISSetIdentity()
263: @*/
264: PetscErrorCode  ISIdentity(IS is,PetscBool  *ident)
265: {

271:   *ident = is->isidentity;
272:   if (*ident) return(0);
273:   if (is->ops->identity) {
274:     (*is->ops->identity)(is,ident);
275:   }
276:   return(0);
277: }

279: /*@
280:    ISSetIdentity - Informs the index set that it is an identity.

282:    Logically Collective on IS

284:    Input Parmeters:
285: .  is - the index set

287:    Level: intermediate


290: .seealso: ISIdentity()
291: @*/
292: PetscErrorCode  ISSetIdentity(IS is)
293: {

298:   is->isidentity = PETSC_TRUE;
299:   ISSetPermutation(is);
300:   return(0);
301: }

303: /*@
304:    ISContiguousLocal - Locates an index set with contiguous range within a global range, if possible

306:    Not Collective

308:    Input Parmeters:
309: +  is - the index set
310: .  gstart - global start
311: -  gend - global end

313:    Output Parameters:
314: +  start - start of contiguous block, as an offset from gstart
315: -  contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE

317:    Level: developer

319: .seealso: ISGetLocalSize(), VecGetOwnershipRange()
320: @*/
321: PetscErrorCode  ISContiguousLocal(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
322: {

329:   if (is->ops->contiguous) {
330:     (*is->ops->contiguous)(is,gstart,gend,start,contig);
331:   } else {
332:     *start  = -1;
333:     *contig = PETSC_FALSE;
334:   }
335:   return(0);
336: }

338: /*@
339:    ISPermutation - PETSC_TRUE or PETSC_FALSE depending on whether the
340:    index set has been declared to be a permutation.

342:    Logically Collective on IS

344:    Input Parmeters:
345: .  is - the index set

347:    Output Parameters:
348: .  perm - PETSC_TRUE if a permutation, else PETSC_FALSE

350:    Level: intermediate


353: .seealso: ISSetPermutation()
354: @*/
355: PetscErrorCode  ISPermutation(IS is,PetscBool  *perm)
356: {
360:   *perm = (PetscBool) is->isperm;
361:   return(0);
362: }

364: /*@
365:    ISSetPermutation - Informs the index set that it is a permutation.

367:    Logically Collective on IS

369:    Input Parmeters:
370: .  is - the index set

372:    Level: intermediate


375:    The debug version of the libraries (./configure --with-debugging=1) checks if the
376:   index set is actually a permutation. The optimized version just believes you.

378: .seealso: ISPermutation()
379: @*/
380: PetscErrorCode  ISSetPermutation(IS is)
381: {
384: #if defined(PETSC_USE_DEBUG)
385:   {
386:     PetscMPIInt    size;

389:     MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
390:     if (size == 1) {
391:       PetscInt       i,n,*idx;
392:       const PetscInt *iidx;

394:       ISGetSize(is,&n);
395:       PetscMalloc1(n,&idx);
396:       ISGetIndices(is,&iidx);
397:       PetscArraycpy(idx,iidx,n);
398:       PetscSortInt(n,idx);
399:       for (i=0; i<n; i++) {
400:         if (idx[i] != i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Index set is not a permutation");
401:       }
402:       PetscFree(idx);
403:       ISRestoreIndices(is,&iidx);
404:     }
405:   }
406: #endif
407:   is->isperm = PETSC_TRUE;
408:   return(0);
409: }

411: /*@
412:    ISDestroy - Destroys an index set.

414:    Collective on IS

416:    Input Parameters:
417: .  is - the index set

419:    Level: beginner

421: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlocked()
422: @*/
423: PetscErrorCode  ISDestroy(IS *is)
424: {

428:   if (!*is) return(0);
430:   if (--((PetscObject)(*is))->refct > 0) {*is = 0; return(0);}
431:   if ((*is)->complement) {
432:     PetscInt refcnt;
433:     PetscObjectGetReference((PetscObject)((*is)->complement), &refcnt);
434:     if (refcnt > 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Nonlocal IS has not been restored");
435:     ISDestroy(&(*is)->complement);
436:   }
437:   if ((*is)->ops->destroy) {
438:     (*(*is)->ops->destroy)(*is);
439:   }
440:   PetscLayoutDestroy(&(*is)->map);
441:   /* Destroy local representations of offproc data. */
442:   PetscFree((*is)->total);
443:   PetscFree((*is)->nonlocal);
444:   PetscHeaderDestroy(is);
445:   return(0);
446: }

448: /*@
449:    ISInvertPermutation - Creates a new permutation that is the inverse of
450:                          a given permutation.

452:    Collective on IS

454:    Input Parameter:
455: +  is - the index set
456: -  nlocal - number of indices on this processor in result (ignored for 1 proccessor) or
457:             use PETSC_DECIDE

459:    Output Parameter:
460: .  isout - the inverse permutation

462:    Level: intermediate

464:    Notes:
465:     For parallel index sets this does the complete parallel permutation, but the
466:     code is not efficient for huge index sets (10,000,000 indices).

468: @*/
469: PetscErrorCode  ISInvertPermutation(IS is,PetscInt nlocal,IS *isout)
470: {

476:   if (!is->isperm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a permutation, must call ISSetPermutation() on the IS first");
477:   if (is->isidentity) {
478:     ISDuplicate(is,isout);
479:   } else {
480:     (*is->ops->invertpermutation)(is,nlocal,isout);
481:     ISSetPermutation(*isout);
482:   }
483:   return(0);
484: }

486: /*@
487:    ISGetSize - Returns the global length of an index set.

489:    Not Collective

491:    Input Parameter:
492: .  is - the index set

494:    Output Parameter:
495: .  size - the global size

497:    Level: beginner


500: @*/
501: PetscErrorCode  ISGetSize(IS is,PetscInt *size)
502: {
506:   *size = is->map->N;
507:   return(0);
508: }

510: /*@
511:    ISGetLocalSize - Returns the local (processor) length of an index set.

513:    Not Collective

515:    Input Parameter:
516: .  is - the index set

518:    Output Parameter:
519: .  size - the local size

521:    Level: beginner

523: @*/
524: PetscErrorCode  ISGetLocalSize(IS is,PetscInt *size)
525: {
529:   *size = is->map->n;
530:   return(0);
531: }

533: /*@C
534:    ISGetIndices - Returns a pointer to the indices.  The user should call
535:    ISRestoreIndices() after having looked at the indices.  The user should
536:    NOT change the indices.

538:    Not Collective

540:    Input Parameter:
541: .  is - the index set

543:    Output Parameter:
544: .  ptr - the location to put the pointer to the indices

546:    Fortran Note:
547:    This routine has two different interfaces from Fortran; the first is not recommend, it does not require Fortran 90
548: $    IS          is
549: $    integer     is_array(1)
550: $    PetscOffset i_is
551: $    int         ierr
552: $       call ISGetIndices(is,is_array,i_is,ierr)
553: $
554: $   Access first local entry in list
555: $      value = is_array(i_is + 1)
556: $
557: $      ...... other code
558: $       call ISRestoreIndices(is,is_array,i_is,ierr)
559:    The second Fortran interface is recommended.
560: $          use petscisdef
561: $          PetscInt, pointer :: array(:)
562: $          PetscErrorCode  ierr
563: $          IS       i
564: $          call ISGetIndicesF90(i,array,ierr)



568:    See the Fortran chapter of the users manual and
569:    petsc/src/is/examples/[tutorials,tests] for details.

571:    Level: intermediate


574: .seealso: ISRestoreIndices(), ISGetIndicesF90()
575: @*/
576: PetscErrorCode  ISGetIndices(IS is,const PetscInt *ptr[])
577: {

583:   (*is->ops->getindices)(is,ptr);
584:   return(0);
585: }

587: /*@C
588:    ISGetMinMax - Gets the minimum and maximum values in an IS

590:    Not Collective

592:    Input Parameter:
593: .  is - the index set

595:    Output Parameter:
596: +   min - the minimum value
597: -   max - the maximum value

599:    Level: intermediate

601:    Notes:
602:     Empty index sets return min=PETSC_MAX_INT and max=PETSC_MIN_INT.
603:     In parallel, it returns the min and max of the local portion of the IS


606: .seealso: ISGetIndices(), ISRestoreIndices(), ISGetIndicesF90()
607: @*/
608: PetscErrorCode  ISGetMinMax(IS is,PetscInt *min,PetscInt *max)
609: {
612:   if (min) *min = is->min;
613:   if (max) *max = is->max;
614:   return(0);
615: }

617: /*@
618:   ISLocate - determine the location of an index within the local component of an index set

620:   Not Collective

622:   Input Parameter:
623: + is - the index set
624: - key - the search key

626:   Output Parameter:
627: . location - if >= 0, a location within the index set that is equal to the key, otherwise the key is not in the index set

629:   Level: intermediate
630: @*/
631: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
632: {

636:   if (is->ops->locate) {
637:     (*is->ops->locate)(is,key,location);
638:   } else {
639:     PetscInt       numIdx;
640:     PetscBool      sorted;
641:     const PetscInt *idx;

643:     ISGetLocalSize(is,&numIdx);
644:     ISGetIndices(is,&idx);
645:     ISSorted(is,&sorted);
646:     if (sorted) {
647:       PetscFindInt(key,numIdx,idx,location);
648:     } else {
649:       PetscInt i;

651:       *location = -1;
652:       for (i = 0; i < numIdx; i++) {
653:         if (idx[i] == key) {
654:           *location = i;
655:           break;
656:         }
657:       }
658:     }
659:     ISRestoreIndices(is,&idx);
660:   }
661:   return(0);
662: }

664: /*@C
665:    ISRestoreIndices - Restores an index set to a usable state after a call
666:                       to ISGetIndices().

668:    Not Collective

670:    Input Parameters:
671: +  is - the index set
672: -  ptr - the pointer obtained by ISGetIndices()

674:    Fortran Note:
675:    This routine is used differently from Fortran
676: $    IS          is
677: $    integer     is_array(1)
678: $    PetscOffset i_is
679: $    int         ierr
680: $       call ISGetIndices(is,is_array,i_is,ierr)
681: $
682: $   Access first local entry in list
683: $      value = is_array(i_is + 1)
684: $
685: $      ...... other code
686: $       call ISRestoreIndices(is,is_array,i_is,ierr)

688:    See the Fortran chapter of the users manual and
689:    petsc/src/is/examples/[tutorials,tests] for details.

691:    Level: intermediate

693:    Note:
694:    This routine zeros out ptr. This is to prevent accidental us of the array after it has been restored.

696: .seealso: ISGetIndices(), ISRestoreIndicesF90()
697: @*/
698: PetscErrorCode  ISRestoreIndices(IS is,const PetscInt *ptr[])
699: {

705:   if (is->ops->restoreindices) {
706:     (*is->ops->restoreindices)(is,ptr);
707:   }
708:   return(0);
709: }

711: static PetscErrorCode ISGatherTotal_Private(IS is)
712: {
714:   PetscInt       i,n,N;
715:   const PetscInt *lindices;
716:   MPI_Comm       comm;
717:   PetscMPIInt    rank,size,*sizes = NULL,*offsets = NULL,nn;


722:   PetscObjectGetComm((PetscObject)is,&comm);
723:   MPI_Comm_size(comm,&size);
724:   MPI_Comm_rank(comm,&rank);
725:   ISGetLocalSize(is,&n);
726:   PetscMalloc2(size,&sizes,size,&offsets);

728:   PetscMPIIntCast(n,&nn);
729:   MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
730:   offsets[0] = 0;
731:   for (i=1; i<size; ++i) offsets[i] = offsets[i-1] + sizes[i-1];
732:   N = offsets[size-1] + sizes[size-1];

734:   PetscMalloc1(N,&(is->total));
735:   ISGetIndices(is,&lindices);
736:   MPI_Allgatherv((void*)lindices,nn,MPIU_INT,is->total,sizes,offsets,MPIU_INT,comm);
737:   ISRestoreIndices(is,&lindices);
738:   is->local_offset = offsets[rank];
739:   PetscFree2(sizes,offsets);
740:   return(0);
741: }

743: /*@C
744:    ISGetTotalIndices - Retrieve an array containing all indices across the communicator.

746:    Collective on IS

748:    Input Parameter:
749: .  is - the index set

751:    Output Parameter:
752: .  indices - total indices with rank 0 indices first, and so on; total array size is
753:              the same as returned with ISGetSize().

755:    Level: intermediate

757:    Notes:
758:     this is potentially nonscalable, but depends on the size of the total index set
759:      and the size of the communicator. This may be feasible for index sets defined on
760:      subcommunicators, such that the set size does not grow with PETSC_WORLD_COMM.
761:      Note also that there is no way to tell where the local part of the indices starts
762:      (use ISGetIndices() and ISGetNonlocalIndices() to retrieve just the local and just
763:       the nonlocal part (complement), respectively).

765: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices(), ISGetSize()
766: @*/
767: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
768: {
770:   PetscMPIInt    size;

775:   MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
776:   if (size == 1) {
777:     (*is->ops->getindices)(is,indices);
778:   } else {
779:     if (!is->total) {
780:       ISGatherTotal_Private(is);
781:     }
782:     *indices = is->total;
783:   }
784:   return(0);
785: }

787: /*@C
788:    ISRestoreTotalIndices - Restore the index array obtained with ISGetTotalIndices().

790:    Not Collective.

792:    Input Parameter:
793: +  is - the index set
794: -  indices - index array; must be the array obtained with ISGetTotalIndices()

796:    Level: intermediate

798: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices()
799: @*/
800: PetscErrorCode  ISRestoreTotalIndices(IS is, const PetscInt *indices[])
801: {
803:   PetscMPIInt    size;

808:   MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
809:   if (size == 1) {
810:     (*is->ops->restoreindices)(is,indices);
811:   } else {
812:     if (is->total != *indices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index array pointer being restored does not point to the array obtained from the IS.");
813:   }
814:   return(0);
815: }
816: /*@C
817:    ISGetNonlocalIndices - Retrieve an array of indices from remote processors
818:                        in this communicator.

820:    Collective on IS

822:    Input Parameter:
823: .  is - the index set

825:    Output Parameter:
826: .  indices - indices with rank 0 indices first, and so on,  omitting
827:              the current rank.  Total number of indices is the difference
828:              total and local, obtained with ISGetSize() and ISGetLocalSize(),
829:              respectively.

831:    Level: intermediate

833:    Notes:
834:     restore the indices using ISRestoreNonlocalIndices().
835:           The same scalability considerations as those for ISGetTotalIndices
836:           apply here.

838: .seealso: ISGetTotalIndices(), ISRestoreNonlocalIndices(), ISGetSize(), ISGetLocalSize().
839: @*/
840: PetscErrorCode  ISGetNonlocalIndices(IS is, const PetscInt *indices[])
841: {
843:   PetscMPIInt    size;
844:   PetscInt       n, N;

849:   MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
850:   if (size == 1) *indices = NULL;
851:   else {
852:     if (!is->total) {
853:       ISGatherTotal_Private(is);
854:     }
855:     ISGetLocalSize(is,&n);
856:     ISGetSize(is,&N);
857:     PetscMalloc1(N-n, &(is->nonlocal));
858:     PetscArraycpy(is->nonlocal, is->total, is->local_offset);
859:     PetscArraycpy(is->nonlocal+is->local_offset, is->total+is->local_offset+n,N - is->local_offset - n);
860:     *indices = is->nonlocal;
861:   }
862:   return(0);
863: }

865: /*@C
866:    ISRestoreTotalIndices - Restore the index array obtained with ISGetNonlocalIndices().

868:    Not Collective.

870:    Input Parameter:
871: +  is - the index set
872: -  indices - index array; must be the array obtained with ISGetNonlocalIndices()

874:    Level: intermediate

876: .seealso: ISGetTotalIndices(), ISGetNonlocalIndices(), ISRestoreTotalIndices()
877: @*/
878: PetscErrorCode  ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
879: {
883:   if (is->nonlocal != *indices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index array pointer being restored does not point to the array obtained from the IS.");
884:   return(0);
885: }

887: /*@
888:    ISGetNonlocalIS - Gather all nonlocal indices for this IS and present
889:                      them as another sequential index set.


892:    Collective on IS

894:    Input Parameter:
895: .  is - the index set

897:    Output Parameter:
898: .  complement - sequential IS with indices identical to the result of
899:                 ISGetNonlocalIndices()

901:    Level: intermediate

903:    Notes:
904:     complement represents the result of ISGetNonlocalIndices as an IS.
905:           Therefore scalability issues similar to ISGetNonlocalIndices apply.
906:           The resulting IS must be restored using ISRestoreNonlocalIS().

908: .seealso: ISGetNonlocalIndices(), ISRestoreNonlocalIndices(),  ISAllGather(), ISGetSize()
909: @*/
910: PetscErrorCode  ISGetNonlocalIS(IS is, IS *complement)
911: {

917:   /* Check if the complement exists already. */
918:   if (is->complement) {
919:     *complement = is->complement;
920:     PetscObjectReference((PetscObject)(is->complement));
921:   } else {
922:     PetscInt       N, n;
923:     const PetscInt *idx;
924:     ISGetSize(is, &N);
925:     ISGetLocalSize(is,&n);
926:     ISGetNonlocalIndices(is, &idx);
927:     ISCreateGeneral(PETSC_COMM_SELF, N-n,idx, PETSC_USE_POINTER, &(is->complement));
928:     PetscObjectReference((PetscObject)is->complement);
929:     *complement = is->complement;
930:   }
931:   return(0);
932: }


935: /*@
936:    ISRestoreNonlocalIS - Restore the IS obtained with ISGetNonlocalIS().

938:    Not collective.

940:    Input Parameter:
941: +  is         - the index set
942: -  complement - index set of is's nonlocal indices

944:    Level: intermediate


947: .seealso: ISGetNonlocalIS(), ISGetNonlocalIndices(), ISRestoreNonlocalIndices()
948: @*/
949: PetscErrorCode  ISRestoreNonlocalIS(IS is, IS *complement)
950: {
952:   PetscInt       refcnt;

957:   if (*complement != is->complement) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Complement IS being restored was not obtained with ISGetNonlocalIS()");
958:   PetscObjectGetReference((PetscObject)(is->complement), &refcnt);
959:   if (refcnt <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Duplicate call to ISRestoreNonlocalIS() detected");
960:   PetscObjectDereference((PetscObject)(is->complement));
961:   return(0);
962: }

964: /*@C
965:    ISView - Displays an index set.

967:    Collective on IS

969:    Input Parameters:
970: +  is - the index set
971: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

973:    Level: intermediate

975: .seealso: PetscViewerASCIIOpen()
976: @*/
977: PetscErrorCode  ISView(IS is,PetscViewer viewer)
978: {

983:   if (!viewer) {
984:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is),&viewer);
985:   }

989:   PetscObjectPrintClassNamePrefixType((PetscObject)is,viewer);
990:   (*is->ops->view)(is,viewer);
991:   return(0);
992: }

994: /*@
995:   ISLoad - Loads a vector that has been stored in binary or HDF5 format with ISView().

997:   Collective on PetscViewer

999:   Input Parameters:
1000: + is - the newly loaded vector, this needs to have been created with ISCreate() or some related function before a call to ISLoad().
1001: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or HDF5 file viewer, obtained from PetscViewerHDF5Open()

1003:   Level: intermediate

1005:   Notes:
1006:   IF using HDF5, you must assign the IS the same name as was used in the IS
1007:   that was stored in the file using PetscObjectSetName(). Otherwise you will
1008:   get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"

1010: .seealso: PetscViewerBinaryOpen(), ISView(), MatLoad(), VecLoad()
1011: @*/
1012: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
1013: {
1014:   PetscBool      isbinary, ishdf5;

1020:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
1021:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
1022:   if (!isbinary && !ishdf5) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1023:   if (!((PetscObject)is)->type_name) {ISSetType(is, ISGENERAL);}
1024:   (*is->ops->load)(is, viewer);
1025:   return(0);
1026: }

1028: /*@
1029:    ISSort - Sorts the indices of an index set.

1031:    Collective on IS

1033:    Input Parameters:
1034: .  is - the index set

1036:    Level: intermediate


1039: .seealso: ISSortRemoveDups(), ISSorted()
1040: @*/
1041: PetscErrorCode  ISSort(IS is)
1042: {

1047:   (*is->ops->sort)(is);
1048:   return(0);
1049: }

1051: /*@
1052:   ISSortRemoveDups - Sorts the indices of an index set, removing duplicates.

1054:   Collective on IS

1056:   Input Parameters:
1057: . is - the index set

1059:   Level: intermediate


1062: .seealso: ISSort(), ISSorted()
1063: @*/
1064: PetscErrorCode ISSortRemoveDups(IS is)
1065: {

1070:   (*is->ops->sortremovedups)(is);
1071:   return(0);
1072: }

1074: /*@
1075:    ISToGeneral - Converts an IS object of any type to ISGENERAL type

1077:    Collective on IS

1079:    Input Parameters:
1080: .  is - the index set

1082:    Level: intermediate


1085: .seealso: ISSorted()
1086: @*/
1087: PetscErrorCode  ISToGeneral(IS is)
1088: {

1093:   if (is->ops->togeneral) {
1094:     (*is->ops->togeneral)(is);
1095:   } else SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Not written for this type %s",((PetscObject)is)->type_name);
1096:   return(0);
1097: }

1099: /*@
1100:    ISSorted - Checks the indices to determine whether they have been sorted.

1102:    Collective on IS

1104:    Input Parameter:
1105: .  is - the index set

1107:    Output Parameter:
1108: .  flg - output flag, either PETSC_TRUE if the index set is sorted,
1109:          or PETSC_FALSE otherwise.

1111:    Notes:
1112:     For parallel IS objects this only indicates if the local part of the IS
1113:           is sorted. So some processors may return PETSC_TRUE while others may
1114:           return PETSC_FALSE.

1116:    Level: intermediate

1118: .seealso: ISSort(), ISSortRemoveDups()
1119: @*/
1120: PetscErrorCode  ISSorted(IS is,PetscBool  *flg)
1121: {

1127:   (*is->ops->sorted)(is,flg);
1128:   return(0);
1129: }

1131: /*@
1132:    ISDuplicate - Creates a duplicate copy of an index set.

1134:    Collective on IS

1136:    Input Parmeters:
1137: .  is - the index set

1139:    Output Parameters:
1140: .  isnew - the copy of the index set

1142:    Level: beginner

1144: .seealso: ISCreateGeneral(), ISCopy()
1145: @*/
1146: PetscErrorCode  ISDuplicate(IS is,IS *newIS)
1147: {

1153:   (*is->ops->duplicate)(is,newIS);
1154:   (*newIS)->isidentity = is->isidentity;
1155:   (*newIS)->isperm     = is->isperm;
1156:   return(0);
1157: }

1159: /*@
1160:    ISCopy - Copies an index set.

1162:    Collective on IS

1164:    Input Parmeters:
1165: .  is - the index set

1167:    Output Parameters:
1168: .  isy - the copy of the index set

1170:    Level: beginner

1172: .seealso: ISDuplicate()
1173: @*/
1174: PetscErrorCode  ISCopy(IS is,IS isy)
1175: {

1182:   if (is == isy) return(0);
1183:   (*is->ops->copy)(is,isy);
1184:   isy->isperm     = is->isperm;
1185:   isy->max        = is->max;
1186:   isy->min        = is->min;
1187:   isy->isidentity = is->isidentity;
1188:   return(0);
1189: }

1191: /*@
1192:    ISOnComm - Split a parallel IS on subcomms (usually self) or concatenate index sets on subcomms into a parallel index set

1194:    Collective on IS

1196:    Input Arguments:
1197: + is - index set
1198: . comm - communicator for new index set
1199: - mode - copy semantics, PETSC_USE_POINTER for no-copy if possible, otherwise PETSC_COPY_VALUES

1201:    Output Arguments:
1202: . newis - new IS on comm

1204:    Level: advanced

1206:    Notes:
1207:    It is usually desirable to create a parallel IS and look at the local part when necessary.

1209:    This function is useful if serial ISs must be created independently, or to view many
1210:    logically independent serial ISs.

1212:    The input IS must have the same type on every process.

1214: .seealso: ISSplit()
1215: @*/
1216: PetscErrorCode  ISOnComm(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
1217: {
1219:   PetscMPIInt    match;

1224:   MPI_Comm_compare(PetscObjectComm((PetscObject)is),comm,&match);
1225:   if (mode != PETSC_COPY_VALUES && (match == MPI_IDENT || match == MPI_CONGRUENT)) {
1226:     PetscObjectReference((PetscObject)is);
1227:     *newis = is;
1228:   } else {
1229:     (*is->ops->oncomm)(is,comm,mode,newis);
1230:   }
1231:   return(0);
1232: }

1234: /*@
1235:    ISSetBlockSize - informs an index set that it has a given block size

1237:    Logicall Collective on IS

1239:    Input Arguments:
1240: + is - index set
1241: - bs - block size

1243:    Level: intermediate

1245: .seealso: ISGetBlockSize(), ISCreateBlock()
1246: @*/
1247: PetscErrorCode  ISSetBlockSize(IS is,PetscInt bs)
1248: {

1254:   if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_OUTOFRANGE,"Block size %D, must be positive",bs);
1255:   (*is->ops->setblocksize)(is,bs);
1256:   return(0);
1257: }

1259: /*@
1260:    ISGetBlockSize - Returns the number of elements in a block.

1262:    Not Collective

1264:    Input Parameter:
1265: .  is - the index set

1267:    Output Parameter:
1268: .  size - the number of elements in a block

1270:    Level: intermediate


1273: .seealso: ISBlockGetSize(), ISGetSize(), ISCreateBlock(), ISSetBlockSize()
1274: @*/
1275: PetscErrorCode  ISGetBlockSize(IS is,PetscInt *size)
1276: {

1280:   PetscLayoutGetBlockSize(is->map, size);
1281:   return(0);
1282: }

1284: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[])
1285: {
1287:   PetscInt       len,i;
1288:   const PetscInt *ptr;

1291:   ISGetSize(is,&len);
1292:   ISGetIndices(is,&ptr);
1293:   for (i=0; i<len; i++) idx[i] = ptr[i];
1294:   ISRestoreIndices(is,&ptr);
1295:   return(0);
1296: }

1298: /*MC
1299:     ISGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1300:     The users should call ISRestoreIndicesF90() after having looked at the
1301:     indices.  The user should NOT change the indices.

1303:     Synopsis:
1304:     ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1306:     Not collective

1308:     Input Parameter:
1309: .   x - index set

1311:     Output Parameters:
1312: +   xx_v - the Fortran90 pointer to the array
1313: -   ierr - error code

1315:     Example of Usage:
1316: .vb
1317:     PetscInt, pointer xx_v(:)
1318:     ....
1319:     call ISGetIndicesF90(x,xx_v,ierr)
1320:     a = xx_v(3)
1321:     call ISRestoreIndicesF90(x,xx_v,ierr)
1322: .ve

1324:     Level: intermediate

1326: .seealso:  ISRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices()


1329: M*/

1331: /*MC
1332:     ISRestoreIndicesF90 - Restores an index set to a usable state after
1333:     a call to ISGetIndicesF90().

1335:     Synopsis:
1336:     ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1338:     Not collective

1340:     Input Parameters:
1341: +   x - index set
1342: -   xx_v - the Fortran90 pointer to the array

1344:     Output Parameter:
1345: .   ierr - error code


1348:     Example of Usage:
1349: .vb
1350:     PetscInt, pointer xx_v(:)
1351:     ....
1352:     call ISGetIndicesF90(x,xx_v,ierr)
1353:     a = xx_v(3)
1354:     call ISRestoreIndicesF90(x,xx_v,ierr)
1355: .ve

1357:     Level: intermediate

1359: .seealso:  ISGetIndicesF90(), ISGetIndices(), ISRestoreIndices()

1361: M*/

1363: /*MC
1364:     ISBlockGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1365:     The users should call ISBlockRestoreIndicesF90() after having looked at the
1366:     indices.  The user should NOT change the indices.

1368:     Synopsis:
1369:     ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1371:     Not collective

1373:     Input Parameter:
1374: .   x - index set

1376:     Output Parameters:
1377: +   xx_v - the Fortran90 pointer to the array
1378: -   ierr - error code
1379:     Example of Usage:
1380: .vb
1381:     PetscInt, pointer xx_v(:)
1382:     ....
1383:     call ISBlockGetIndicesF90(x,xx_v,ierr)
1384:     a = xx_v(3)
1385:     call ISBlockRestoreIndicesF90(x,xx_v,ierr)
1386: .ve

1388:     Level: intermediate

1390: .seealso:  ISBlockRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices(),
1391:            ISRestoreIndices()

1393: M*/

1395: /*MC
1396:     ISBlockRestoreIndicesF90 - Restores an index set to a usable state after
1397:     a call to ISBlockGetIndicesF90().

1399:     Synopsis:
1400:     ISBlockRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1402:     Not Collective

1404:     Input Parameters:
1405: +   x - index set
1406: -   xx_v - the Fortran90 pointer to the array

1408:     Output Parameter:
1409: .   ierr - error code

1411:     Example of Usage:
1412: .vb
1413:     PetscInt, pointer xx_v(:)
1414:     ....
1415:     call ISBlockGetIndicesF90(x,xx_v,ierr)
1416:     a = xx_v(3)
1417:     call ISBlockRestoreIndicesF90(x,xx_v,ierr)
1418: .ve

1420:     Notes:
1421:     Not yet supported for all F90 compilers

1423:     Level: intermediate

1425: .seealso:  ISBlockGetIndicesF90(), ISGetIndices(), ISRestoreIndices(), ISRestoreIndicesF90()

1427: M*/