Actual source code: index.c

petsc-master 2019-07-18
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: {

508:   (*is->ops->getsize)(is,size);
509:   return(0);
510: }

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

515:    Not Collective

517:    Input Parameter:
518: .  is - the index set

520:    Output Parameter:
521: .  size - the local size

523:    Level: beginner

525: @*/
526: PetscErrorCode  ISGetLocalSize(IS is,PetscInt *size)
527: {

533:   (*is->ops->getlocalsize)(is,size);
534:   return(0);
535: }

537: /*@C
538:    ISGetIndices - Returns a pointer to the indices.  The user should call
539:    ISRestoreIndices() after having looked at the indices.  The user should
540:    NOT change the indices.

542:    Not Collective

544:    Input Parameter:
545: .  is - the index set

547:    Output Parameter:
548: .  ptr - the location to put the pointer to the indices

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



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

575:    Level: intermediate


578: .seealso: ISRestoreIndices(), ISGetIndicesF90()
579: @*/
580: PetscErrorCode  ISGetIndices(IS is,const PetscInt *ptr[])
581: {

587:   (*is->ops->getindices)(is,ptr);
588:   return(0);
589: }

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

594:    Not Collective

596:    Input Parameter:
597: .  is - the index set

599:    Output Parameter:
600: +   min - the minimum value
601: -   max - the maximum value

603:    Level: intermediate

605:    Notes:
606:     Empty index sets return min=PETSC_MAX_INT and max=PETSC_MIN_INT.
607:     In parallel, it returns the min and max of the local portion of the IS


610: .seealso: ISGetIndices(), ISRestoreIndices(), ISGetIndicesF90()
611: @*/
612: PetscErrorCode  ISGetMinMax(IS is,PetscInt *min,PetscInt *max)
613: {
616:   if (min) *min = is->min;
617:   if (max) *max = is->max;
618:   return(0);
619: }

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

624:   Not Collective

626:   Input Parameter:
627: + is - the index set
628: - key - the search key

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

633:   Level: intermediate
634: @*/
635: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
636: {

640:   if (is->ops->locate) {
641:     (*is->ops->locate)(is,key,location);
642:   } else {
643:     PetscInt       numIdx;
644:     PetscBool      sorted;
645:     const PetscInt *idx;

647:     ISGetLocalSize(is,&numIdx);
648:     ISGetIndices(is,&idx);
649:     ISSorted(is,&sorted);
650:     if (sorted) {
651:       PetscFindInt(key,numIdx,idx,location);
652:     } else {
653:       PetscInt i;

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

668: /*@C
669:    ISRestoreIndices - Restores an index set to a usable state after a call
670:                       to ISGetIndices().

672:    Not Collective

674:    Input Parameters:
675: +  is - the index set
676: -  ptr - the pointer obtained by ISGetIndices()

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

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

695:    Level: intermediate

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

700: .seealso: ISGetIndices(), ISRestoreIndicesF90()
701: @*/
702: PetscErrorCode  ISRestoreIndices(IS is,const PetscInt *ptr[])
703: {

709:   if (is->ops->restoreindices) {
710:     (*is->ops->restoreindices)(is,ptr);
711:   }
712:   return(0);
713: }

715: static PetscErrorCode ISGatherTotal_Private(IS is)
716: {
718:   PetscInt       i,n,N;
719:   const PetscInt *lindices;
720:   MPI_Comm       comm;
721:   PetscMPIInt    rank,size,*sizes = NULL,*offsets = NULL,nn;


726:   PetscObjectGetComm((PetscObject)is,&comm);
727:   MPI_Comm_size(comm,&size);
728:   MPI_Comm_rank(comm,&rank);
729:   ISGetLocalSize(is,&n);
730:   PetscMalloc2(size,&sizes,size,&offsets);

732:   PetscMPIIntCast(n,&nn);
733:   MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
734:   offsets[0] = 0;
735:   for (i=1; i<size; ++i) offsets[i] = offsets[i-1] + sizes[i-1];
736:   N = offsets[size-1] + sizes[size-1];

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

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

750:    Collective on IS

752:    Input Parameter:
753: .  is - the index set

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

759:    Level: intermediate

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

769: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices(), ISGetSize()
770: @*/
771: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
772: {
774:   PetscMPIInt    size;

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

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

794:    Not Collective.

796:    Input Parameter:
797: +  is - the index set
798: -  indices - index array; must be the array obtained with ISGetTotalIndices()

800:    Level: intermediate

802: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices()
803: @*/
804: PetscErrorCode  ISRestoreTotalIndices(IS is, const PetscInt *indices[])
805: {
807:   PetscMPIInt    size;

812:   MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
813:   if (size == 1) {
814:     (*is->ops->restoreindices)(is,indices);
815:   } else {
816:     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.");
817:   }
818:   return(0);
819: }
820: /*@C
821:    ISGetNonlocalIndices - Retrieve an array of indices from remote processors
822:                        in this communicator.

824:    Collective on IS

826:    Input Parameter:
827: .  is - the index set

829:    Output Parameter:
830: .  indices - indices with rank 0 indices first, and so on,  omitting
831:              the current rank.  Total number of indices is the difference
832:              total and local, obtained with ISGetSize() and ISGetLocalSize(),
833:              respectively.

835:    Level: intermediate

837:    Notes:
838:     restore the indices using ISRestoreNonlocalIndices().
839:           The same scalability considerations as those for ISGetTotalIndices
840:           apply here.

842: .seealso: ISGetTotalIndices(), ISRestoreNonlocalIndices(), ISGetSize(), ISGetLocalSize().
843: @*/
844: PetscErrorCode  ISGetNonlocalIndices(IS is, const PetscInt *indices[])
845: {
847:   PetscMPIInt    size;
848:   PetscInt       n, N;

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

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

872:    Not Collective.

874:    Input Parameter:
875: +  is - the index set
876: -  indices - index array; must be the array obtained with ISGetNonlocalIndices()

878:    Level: intermediate

880: .seealso: ISGetTotalIndices(), ISGetNonlocalIndices(), ISRestoreTotalIndices()
881: @*/
882: PetscErrorCode  ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
883: {
887:   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.");
888:   return(0);
889: }

891: /*@
892:    ISGetNonlocalIS - Gather all nonlocal indices for this IS and present
893:                      them as another sequential index set.


896:    Collective on IS

898:    Input Parameter:
899: .  is - the index set

901:    Output Parameter:
902: .  complement - sequential IS with indices identical to the result of
903:                 ISGetNonlocalIndices()

905:    Level: intermediate

907:    Notes:
908:     complement represents the result of ISGetNonlocalIndices as an IS.
909:           Therefore scalability issues similar to ISGetNonlocalIndices apply.
910:           The resulting IS must be restored using ISRestoreNonlocalIS().

912: .seealso: ISGetNonlocalIndices(), ISRestoreNonlocalIndices(),  ISAllGather(), ISGetSize()
913: @*/
914: PetscErrorCode  ISGetNonlocalIS(IS is, IS *complement)
915: {

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


939: /*@
940:    ISRestoreNonlocalIS - Restore the IS obtained with ISGetNonlocalIS().

942:    Not collective.

944:    Input Parameter:
945: +  is         - the index set
946: -  complement - index set of is's nonlocal indices

948:    Level: intermediate


951: .seealso: ISGetNonlocalIS(), ISGetNonlocalIndices(), ISRestoreNonlocalIndices()
952: @*/
953: PetscErrorCode  ISRestoreNonlocalIS(IS is, IS *complement)
954: {
956:   PetscInt       refcnt;

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

968: /*@C
969:    ISView - Displays an index set.

971:    Collective on IS

973:    Input Parameters:
974: +  is - the index set
975: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

977:    Level: intermediate

979: .seealso: PetscViewerASCIIOpen()
980: @*/
981: PetscErrorCode  ISView(IS is,PetscViewer viewer)
982: {

987:   if (!viewer) {
988:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is),&viewer);
989:   }

993:   PetscObjectPrintClassNamePrefixType((PetscObject)is,viewer);
994:   (*is->ops->view)(is,viewer);
995:   return(0);
996: }

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

1001:   Collective on PetscViewer

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

1007:   Level: intermediate

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

1014: .seealso: PetscViewerBinaryOpen(), ISView(), MatLoad(), VecLoad()
1015: @*/
1016: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
1017: {
1018:   PetscBool      isbinary, ishdf5;

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

1032: /*@
1033:    ISSort - Sorts the indices of an index set.

1035:    Collective on IS

1037:    Input Parameters:
1038: .  is - the index set

1040:    Level: intermediate


1043: .seealso: ISSortRemoveDups(), ISSorted()
1044: @*/
1045: PetscErrorCode  ISSort(IS is)
1046: {

1051:   (*is->ops->sort)(is);
1052:   return(0);
1053: }

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

1058:   Collective on IS

1060:   Input Parameters:
1061: . is - the index set

1063:   Level: intermediate


1066: .seealso: ISSort(), ISSorted()
1067: @*/
1068: PetscErrorCode ISSortRemoveDups(IS is)
1069: {

1074:   (*is->ops->sortremovedups)(is);
1075:   return(0);
1076: }

1078: /*@
1079:    ISToGeneral - Converts an IS object of any type to ISGENERAL type

1081:    Collective on IS

1083:    Input Parameters:
1084: .  is - the index set

1086:    Level: intermediate


1089: .seealso: ISSorted()
1090: @*/
1091: PetscErrorCode  ISToGeneral(IS is)
1092: {

1097:   if (is->ops->togeneral) {
1098:     (*is->ops->togeneral)(is);
1099:   } else SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Not written for this type %s",((PetscObject)is)->type_name);
1100:   return(0);
1101: }

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

1106:    Collective on IS

1108:    Input Parameter:
1109: .  is - the index set

1111:    Output Parameter:
1112: .  flg - output flag, either PETSC_TRUE if the index set is sorted,
1113:          or PETSC_FALSE otherwise.

1115:    Notes:
1116:     For parallel IS objects this only indicates if the local part of the IS
1117:           is sorted. So some processors may return PETSC_TRUE while others may
1118:           return PETSC_FALSE.

1120:    Level: intermediate

1122: .seealso: ISSort(), ISSortRemoveDups()
1123: @*/
1124: PetscErrorCode  ISSorted(IS is,PetscBool  *flg)
1125: {

1131:   (*is->ops->sorted)(is,flg);
1132:   return(0);
1133: }

1135: /*@
1136:    ISDuplicate - Creates a duplicate copy of an index set.

1138:    Collective on IS

1140:    Input Parmeters:
1141: .  is - the index set

1143:    Output Parameters:
1144: .  isnew - the copy of the index set

1146:    Level: beginner

1148: .seealso: ISCreateGeneral(), ISCopy()
1149: @*/
1150: PetscErrorCode  ISDuplicate(IS is,IS *newIS)
1151: {

1157:   (*is->ops->duplicate)(is,newIS);
1158:   (*newIS)->isidentity = is->isidentity;
1159:   (*newIS)->isperm     = is->isperm;
1160:   return(0);
1161: }

1163: /*@
1164:    ISCopy - Copies an index set.

1166:    Collective on IS

1168:    Input Parmeters:
1169: .  is - the index set

1171:    Output Parameters:
1172: .  isy - the copy of the index set

1174:    Level: beginner

1176: .seealso: ISDuplicate()
1177: @*/
1178: PetscErrorCode  ISCopy(IS is,IS isy)
1179: {

1186:   if (is == isy) return(0);
1187:   (*is->ops->copy)(is,isy);
1188:   isy->isperm     = is->isperm;
1189:   isy->max        = is->max;
1190:   isy->min        = is->min;
1191:   isy->isidentity = is->isidentity;
1192:   return(0);
1193: }

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

1198:    Collective on IS

1200:    Input Arguments:
1201: + is - index set
1202: . comm - communicator for new index set
1203: - mode - copy semantics, PETSC_USE_POINTER for no-copy if possible, otherwise PETSC_COPY_VALUES

1205:    Output Arguments:
1206: . newis - new IS on comm

1208:    Level: advanced

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

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

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

1218: .seealso: ISSplit()
1219: @*/
1220: PetscErrorCode  ISOnComm(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
1221: {
1223:   PetscMPIInt    match;

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

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

1241:    Logicall Collective on IS

1243:    Input Arguments:
1244: + is - index set
1245: - bs - block size

1247:    Level: intermediate

1249: .seealso: ISGetBlockSize(), ISCreateBlock()
1250: @*/
1251: PetscErrorCode  ISSetBlockSize(IS is,PetscInt bs)
1252: {

1258:   if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_OUTOFRANGE,"Block size %D, must be positive",bs);
1259:   (*is->ops->setblocksize)(is,bs);
1260:   return(0);
1261: }

1263: /*@
1264:    ISGetBlockSize - Returns the number of elements in a block.

1266:    Not Collective

1268:    Input Parameter:
1269: .  is - the index set

1271:    Output Parameter:
1272: .  size - the number of elements in a block

1274:    Level: intermediate


1277: .seealso: ISBlockGetSize(), ISGetSize(), ISCreateBlock(), ISSetBlockSize()
1278: @*/
1279: PetscErrorCode  ISGetBlockSize(IS is,PetscInt *size)
1280: {

1284:   PetscLayoutGetBlockSize(is->map, size);
1285:   return(0);
1286: }

1288: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[])
1289: {
1291:   PetscInt       len,i;
1292:   const PetscInt *ptr;

1295:   ISGetSize(is,&len);
1296:   ISGetIndices(is,&ptr);
1297:   for (i=0; i<len; i++) idx[i] = ptr[i];
1298:   ISRestoreIndices(is,&ptr);
1299:   return(0);
1300: }

1302: /*MC
1303:     ISGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1304:     The users should call ISRestoreIndicesF90() after having looked at the
1305:     indices.  The user should NOT change the indices.

1307:     Synopsis:
1308:     ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1310:     Not collective

1312:     Input Parameter:
1313: .   x - index set

1315:     Output Parameters:
1316: +   xx_v - the Fortran90 pointer to the array
1317: -   ierr - error code

1319:     Example of Usage:
1320: .vb
1321:     PetscInt, pointer xx_v(:)
1322:     ....
1323:     call ISGetIndicesF90(x,xx_v,ierr)
1324:     a = xx_v(3)
1325:     call ISRestoreIndicesF90(x,xx_v,ierr)
1326: .ve

1328:     Level: intermediate

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


1333: M*/

1335: /*MC
1336:     ISRestoreIndicesF90 - Restores an index set to a usable state after
1337:     a call to ISGetIndicesF90().

1339:     Synopsis:
1340:     ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1342:     Not collective

1344:     Input Parameters:
1345: +   x - index set
1346: -   xx_v - the Fortran90 pointer to the array

1348:     Output Parameter:
1349: .   ierr - error code


1352:     Example of Usage:
1353: .vb
1354:     PetscInt, pointer xx_v(:)
1355:     ....
1356:     call ISGetIndicesF90(x,xx_v,ierr)
1357:     a = xx_v(3)
1358:     call ISRestoreIndicesF90(x,xx_v,ierr)
1359: .ve

1361:     Level: intermediate

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

1365: M*/

1367: /*MC
1368:     ISBlockGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1369:     The users should call ISBlockRestoreIndicesF90() after having looked at the
1370:     indices.  The user should NOT change the indices.

1372:     Synopsis:
1373:     ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1375:     Not collective

1377:     Input Parameter:
1378: .   x - index set

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

1392:     Level: intermediate

1394: .seealso:  ISBlockRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices(),
1395:            ISRestoreIndices()

1397: M*/

1399: /*MC
1400:     ISBlockRestoreIndicesF90 - Restores an index set to a usable state after
1401:     a call to ISBlockGetIndicesF90().

1403:     Synopsis:
1404:     ISBlockRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1406:     Not Collective

1408:     Input Parameters:
1409: +   x - index set
1410: -   xx_v - the Fortran90 pointer to the array

1412:     Output Parameter:
1413: .   ierr - error code

1415:     Example of Usage:
1416: .vb
1417:     PetscInt, pointer xx_v(:)
1418:     ....
1419:     call ISBlockGetIndicesF90(x,xx_v,ierr)
1420:     a = xx_v(3)
1421:     call ISBlockRestoreIndicesF90(x,xx_v,ierr)
1422: .ve

1424:     Notes:
1425:     Not yet supported for all F90 compilers

1427:     Level: intermediate

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

1431: M*/