Actual source code: index.c

petsc-3.9.4 2018-09-11
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;

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

 65:   /* create sf : leaf_data == multiplicity of indexes, root data == global index in layout */
 66:   PetscMalloc2(n,&leaf_data,Nl,&root_data);
 67:   if (subset_mult) {
 68:     const PetscInt* idxs_mult;

 70:     ISGetIndices(subset_mult,&idxs_mult);
 71:     PetscMemcpy(leaf_data,idxs_mult,n*sizeof(PetscInt));
 72:     ISRestoreIndices(subset_mult,&idxs_mult);
 73:   } else {
 74:     for (i=0;i<n;i++) leaf_data[i] = 1;
 75:   }
 76:   /* local size of new subset */
 77:   n_n = 0;
 78:   for (i=0;i<n;i++) n_n += leaf_data[i];

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

 88:   /* reduce from leaves to roots */
 89:   PetscMemzero(root_data,Nl*sizeof(PetscInt));
 90:   PetscSFReduceBegin(sf,MPIU_INT,leaf_data,root_data,MPI_MAX);
 91:   PetscSFReduceEnd(sf,MPIU_INT,leaf_data,root_data,MPI_MAX);

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

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

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

120:   /* adapt root data with cumulative */
121:   if (first_found) {
122:     PetscInt old_index;

124:     root_data[first_index] += start;
125:     old_index = first_index;
126:     for (i=first_index+1;i<Nl;i++) {
127:       if (root_data[i]) {
128:         root_data[i] += root_data[old_index];
129:         old_index = i;
130:       }
131:     }
132:   }

134:   /* from roots to leaves */
135:   PetscSFBcastBegin(sf,MPIU_INT,root_data,leaf_data);
136:   PetscSFBcastEnd(sf,MPIU_INT,root_data,leaf_data);
137:   PetscSFDestroy(&sf);

139:   /* create new IS with global indexes without holes */
140:   if (subset_mult) {
141:     const PetscInt* idxs_mult;
142:     PetscInt        cum;

144:     cum = 0;
145:     ISGetIndices(subset_mult,&idxs_mult);
146:     for (i=0;i<n;i++) {
147:       PetscInt j;
148:       for (j=0;j<idxs_mult[i];j++) gidxs[cum++] = leaf_data[i] - idxs_mult[i] + j;
149:     }
150:     ISRestoreIndices(subset_mult,&idxs_mult);
151:   } else {
152:     for (i=0;i<n;i++) {
153:       gidxs[i] = leaf_data[i]-1;
154:     }
155:   }
156:   ISCreateGeneral(PetscObjectComm((PetscObject)subset),n_n,gidxs,PETSC_OWN_POINTER,subset_n);
157:   PetscFree2(leaf_data,root_data);
158:   return(0);
159: }

161: /*@
162:    ISIdentity - Determines whether index set is the identity mapping.

164:    Collective on IS

166:    Input Parmeters:
167: .  is - the index set

169:    Output Parameters:
170: .  ident - PETSC_TRUE if an identity, else PETSC_FALSE

172:    Level: intermediate

174:    Concepts: identity mapping
175:    Concepts: index sets^is identity

177: .seealso: ISSetIdentity()
178: @*/
179: PetscErrorCode  ISIdentity(IS is,PetscBool  *ident)
180: {

186:   *ident = is->isidentity;
187:   if (*ident) return(0);
188:   if (is->ops->identity) {
189:     (*is->ops->identity)(is,ident);
190:   }
191:   return(0);
192: }

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

197:    Logically Collective on IS

199:    Input Parmeters:
200: .  is - the index set

202:    Level: intermediate

204:    Concepts: identity mapping
205:    Concepts: index sets^is identity

207: .seealso: ISIdentity()
208: @*/
209: PetscErrorCode  ISSetIdentity(IS is)
210: {

215:   is->isidentity = PETSC_TRUE;
216:   ISSetPermutation(is);
217:   return(0);
218: }

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

223:    Not Collective

225:    Input Parmeters:
226: +  is - the index set
227: .  gstart - global start
228: .  gend - global end

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

234:    Level: developer

236:    Concepts: index sets^is contiguous

238: .seealso: ISGetLocalSize(), VecGetOwnershipRange()
239: @*/
240: PetscErrorCode  ISContiguousLocal(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
241: {

248:   if (is->ops->contiguous) {
249:     (*is->ops->contiguous)(is,gstart,gend,start,contig);
250:   } else {
251:     *start  = -1;
252:     *contig = PETSC_FALSE;
253:   }
254:   return(0);
255: }

257: /*@
258:    ISPermutation - PETSC_TRUE or PETSC_FALSE depending on whether the
259:    index set has been declared to be a permutation.

261:    Logically Collective on IS

263:    Input Parmeters:
264: .  is - the index set

266:    Output Parameters:
267: .  perm - PETSC_TRUE if a permutation, else PETSC_FALSE

269:    Level: intermediate

271:   Concepts: permutation
272:   Concepts: index sets^is permutation

274: .seealso: ISSetPermutation()
275: @*/
276: PetscErrorCode  ISPermutation(IS is,PetscBool  *perm)
277: {
281:   *perm = (PetscBool) is->isperm;
282:   return(0);
283: }

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

288:    Logically Collective on IS

290:    Input Parmeters:
291: .  is - the index set

293:    Level: intermediate

295:   Concepts: permutation
296:   Concepts: index sets^permutation

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

301: .seealso: ISPermutation()
302: @*/
303: PetscErrorCode  ISSetPermutation(IS is)
304: {
307: #if defined(PETSC_USE_DEBUG)
308:   {
309:     PetscMPIInt    size;

312:     MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
313:     if (size == 1) {
314:       PetscInt       i,n,*idx;
315:       const PetscInt *iidx;

317:       ISGetSize(is,&n);
318:       PetscMalloc1(n,&idx);
319:       ISGetIndices(is,&iidx);
320:       PetscMemcpy(idx,iidx,n*sizeof(PetscInt));
321:       PetscSortInt(n,idx);
322:       for (i=0; i<n; i++) {
323:         if (idx[i] != i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Index set is not a permutation");
324:       }
325:       PetscFree(idx);
326:       ISRestoreIndices(is,&iidx);
327:     }
328:   }
329: #endif
330:   is->isperm = PETSC_TRUE;
331:   return(0);
332: }

334: /*@
335:    ISDestroy - Destroys an index set.

337:    Collective on IS

339:    Input Parameters:
340: .  is - the index set

342:    Level: beginner

344: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlocked()
345: @*/
346: PetscErrorCode  ISDestroy(IS *is)
347: {

351:   if (!*is) return(0);
353:   if (--((PetscObject)(*is))->refct > 0) {*is = 0; return(0);}
354:   if ((*is)->complement) {
355:     PetscInt refcnt;
356:     PetscObjectGetReference((PetscObject)((*is)->complement), &refcnt);
357:     if (refcnt > 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Nonlocal IS has not been restored");
358:     ISDestroy(&(*is)->complement);
359:   }
360:   if ((*is)->ops->destroy) {
361:     (*(*is)->ops->destroy)(*is);
362:   }
363:   PetscLayoutDestroy(&(*is)->map);
364:   /* Destroy local representations of offproc data. */
365:   PetscFree((*is)->total);
366:   PetscFree((*is)->nonlocal);
367:   PetscHeaderDestroy(is);
368:   return(0);
369: }

371: /*@
372:    ISInvertPermutation - Creates a new permutation that is the inverse of
373:                          a given permutation.

375:    Collective on IS

377:    Input Parameter:
378: +  is - the index set
379: -  nlocal - number of indices on this processor in result (ignored for 1 proccessor) or
380:             use PETSC_DECIDE

382:    Output Parameter:
383: .  isout - the inverse permutation

385:    Level: intermediate

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

390:    Concepts: inverse permutation
391:    Concepts: permutation^inverse
392:    Concepts: index sets^inverting
393: @*/
394: PetscErrorCode  ISInvertPermutation(IS is,PetscInt nlocal,IS *isout)
395: {

401:   if (!is->isperm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a permutation, must call ISSetPermutation() on the IS first");
402:   if (is->isidentity) {
403:     ISDuplicate(is,isout);
404:   } else {
405:     (*is->ops->invertpermutation)(is,nlocal,isout);
406:     ISSetPermutation(*isout);
407:   }
408:   return(0);
409: }

411: /*@
412:    ISGetSize - Returns the global length of an index set.

414:    Not Collective

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

419:    Output Parameter:
420: .  size - the global size

422:    Level: beginner

424:    Concepts: size^of index set
425:    Concepts: index sets^size

427: @*/
428: PetscErrorCode  ISGetSize(IS is,PetscInt *size)
429: {

435:   (*is->ops->getsize)(is,size);
436:   return(0);
437: }

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

442:    Not Collective

444:    Input Parameter:
445: .  is - the index set

447:    Output Parameter:
448: .  size - the local size

450:    Level: beginner

452:    Concepts: size^of index set
453:    Concepts: local size^of index set
454:    Concepts: index sets^local size

456: @*/
457: PetscErrorCode  ISGetLocalSize(IS is,PetscInt *size)
458: {

464:   (*is->ops->getlocalsize)(is,size);
465:   return(0);
466: }

468: /*@C
469:    ISGetIndices - Returns a pointer to the indices.  The user should call
470:    ISRestoreIndices() after having looked at the indices.  The user should
471:    NOT change the indices.

473:    Not Collective

475:    Input Parameter:
476: .  is - the index set

478:    Output Parameter:
479: .  ptr - the location to put the pointer to the indices

481:    Fortran Note:
482:    This routine has two different interfaces from Fortran; the first is not recommend, it does not require Fortran 90
483: $    IS          is
484: $    integer     is_array(1)
485: $    PetscOffset i_is
486: $    int         ierr
487: $       call ISGetIndices(is,is_array,i_is,ierr)
488: $
489: $   Access first local entry in list
490: $      value = is_array(i_is + 1)
491: $
492: $      ...... other code
493: $       call ISRestoreIndices(is,is_array,i_is,ierr)
494:    The second Fortran interface is recommended.
495: $          use petscisdef
496: $          PetscInt, pointer :: array(:)
497: $          PetscErrorCode  ierr
498: $          IS       i
499: $          call ISGetIndicesF90(i,array,ierr)



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

506:    Level: intermediate

508:    Concepts: index sets^getting indices
509:    Concepts: indices of index set

511: .seealso: ISRestoreIndices(), ISGetIndicesF90()
512: @*/
513: PetscErrorCode  ISGetIndices(IS is,const PetscInt *ptr[])
514: {

520:   (*is->ops->getindices)(is,ptr);
521:   return(0);
522: }

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

527:    Not Collective

529:    Input Parameter:
530: .  is - the index set

532:    Output Parameter:
533: +   min - the minimum value
534: -   max - the maximum value

536:    Level: intermediate

538:    Notes:
539:     Empty index sets return min=PETSC_MAX_INT and max=PETSC_MIN_INT.
540:     In parallel, it returns the min and max of the local portion of the IS

542:    Concepts: index sets^getting indices
543:    Concepts: indices of index set

545: .seealso: ISGetIndices(), ISRestoreIndices(), ISGetIndicesF90()
546: @*/
547: PetscErrorCode  ISGetMinMax(IS is,PetscInt *min,PetscInt *max)
548: {
551:   if (min) *min = is->min;
552:   if (max) *max = is->max;
553:   return(0);
554: }

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

559:   Not Collective

561:   Input Parameter:
562: + is - the index set
563: - key - the search key

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

568:   Level: intermediate
569: @*/
570: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
571: {

575:   if (is->ops->locate) {
576:     (*is->ops->locate)(is,key,location);
577:   } else {
578:     PetscInt       numIdx;
579:     PetscBool      sorted;
580:     const PetscInt *idx;

582:     ISGetLocalSize(is,&numIdx);
583:     ISGetIndices(is,&idx);
584:     ISSorted(is,&sorted);
585:     if (sorted) {
586:       PetscFindInt(key,numIdx,idx,location);
587:     } else {
588:       PetscInt i;

590:       *location = -1;
591:       for (i = 0; i < numIdx; i++) {
592:         if (idx[i] == key) {
593:           *location = i;
594:           break;
595:         }
596:       }
597:     }
598:     ISRestoreIndices(is,&idx);
599:   }
600:   return(0);
601: }

603: /*@C
604:    ISRestoreIndices - Restores an index set to a usable state after a call
605:                       to ISGetIndices().

607:    Not Collective

609:    Input Parameters:
610: +  is - the index set
611: -  ptr - the pointer obtained by ISGetIndices()

613:    Fortran Note:
614:    This routine is used differently from Fortran
615: $    IS          is
616: $    integer     is_array(1)
617: $    PetscOffset i_is
618: $    int         ierr
619: $       call ISGetIndices(is,is_array,i_is,ierr)
620: $
621: $   Access first local entry in list
622: $      value = is_array(i_is + 1)
623: $
624: $      ...... other code
625: $       call ISRestoreIndices(is,is_array,i_is,ierr)

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

630:    Level: intermediate

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

635: .seealso: ISGetIndices(), ISRestoreIndicesF90()
636: @*/
637: PetscErrorCode  ISRestoreIndices(IS is,const PetscInt *ptr[])
638: {

644:   if (is->ops->restoreindices) {
645:     (*is->ops->restoreindices)(is,ptr);
646:   }
647:   return(0);
648: }

650: static PetscErrorCode ISGatherTotal_Private(IS is)
651: {
653:   PetscInt       i,n,N;
654:   const PetscInt *lindices;
655:   MPI_Comm       comm;
656:   PetscMPIInt    rank,size,*sizes = NULL,*offsets = NULL,nn;


661:   PetscObjectGetComm((PetscObject)is,&comm);
662:   MPI_Comm_size(comm,&size);
663:   MPI_Comm_rank(comm,&rank);
664:   ISGetLocalSize(is,&n);
665:   PetscMalloc2(size,&sizes,size,&offsets);

667:   PetscMPIIntCast(n,&nn);
668:   MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
669:   offsets[0] = 0;
670:   for (i=1; i<size; ++i) offsets[i] = offsets[i-1] + sizes[i-1];
671:   N = offsets[size-1] + sizes[size-1];

673:   PetscMalloc1(N,&(is->total));
674:   ISGetIndices(is,&lindices);
675:   MPI_Allgatherv((void*)lindices,nn,MPIU_INT,is->total,sizes,offsets,MPIU_INT,comm);
676:   ISRestoreIndices(is,&lindices);
677:   is->local_offset = offsets[rank];
678:   PetscFree2(sizes,offsets);
679:   return(0);
680: }

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

685:    Collective on IS

687:    Input Parameter:
688: .  is - the index set

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

694:    Level: intermediate

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

703:    Concepts: index sets^getting nonlocal indices
704: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices(), ISGetSize()
705: @*/
706: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
707: {
709:   PetscMPIInt    size;

714:   MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
715:   if (size == 1) {
716:     (*is->ops->getindices)(is,indices);
717:   } else {
718:     if (!is->total) {
719:       ISGatherTotal_Private(is);
720:     }
721:     *indices = is->total;
722:   }
723:   return(0);
724: }

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

729:    Not Collective.

731:    Input Parameter:
732: +  is - the index set
733: -  indices - index array; must be the array obtained with ISGetTotalIndices()

735:    Level: intermediate

737:    Concepts: index sets^getting nonlocal indices
738:    Concepts: index sets^restoring nonlocal indices
739: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices()
740: @*/
741: PetscErrorCode  ISRestoreTotalIndices(IS is, const PetscInt *indices[])
742: {
744:   PetscMPIInt    size;

749:   MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
750:   if (size == 1) {
751:     (*is->ops->restoreindices)(is,indices);
752:   } else {
753:     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.");
754:   }
755:   return(0);
756: }
757: /*@C
758:    ISGetNonlocalIndices - Retrieve an array of indices from remote processors
759:                        in this communicator.

761:    Collective on IS

763:    Input Parameter:
764: .  is - the index set

766:    Output Parameter:
767: .  indices - indices with rank 0 indices first, and so on,  omitting
768:              the current rank.  Total number of indices is the difference
769:              total and local, obtained with ISGetSize() and ISGetLocalSize(),
770:              respectively.

772:    Level: intermediate

774:    Notes: restore the indices using ISRestoreNonlocalIndices().
775:           The same scalability considerations as those for ISGetTotalIndices
776:           apply here.

778:    Concepts: index sets^getting nonlocal indices
779: .seealso: ISGetTotalIndices(), ISRestoreNonlocalIndices(), ISGetSize(), ISGetLocalSize().
780: @*/
781: PetscErrorCode  ISGetNonlocalIndices(IS is, const PetscInt *indices[])
782: {
784:   PetscMPIInt    size;
785:   PetscInt       n, N;

790:   MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
791:   if (size == 1) *indices = NULL;
792:   else {
793:     if (!is->total) {
794:       ISGatherTotal_Private(is);
795:     }
796:     ISGetLocalSize(is,&n);
797:     ISGetSize(is,&N);
798:     PetscMalloc1(N-n, &(is->nonlocal));
799:     PetscMemcpy(is->nonlocal, is->total, sizeof(PetscInt)*is->local_offset);
800:     PetscMemcpy(is->nonlocal+is->local_offset, is->total+is->local_offset+n, sizeof(PetscInt)*(N - is->local_offset - n));
801:     *indices = is->nonlocal;
802:   }
803:   return(0);
804: }

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

809:    Not Collective.

811:    Input Parameter:
812: +  is - the index set
813: -  indices - index array; must be the array obtained with ISGetNonlocalIndices()

815:    Level: intermediate

817:    Concepts: index sets^getting nonlocal indices
818:    Concepts: index sets^restoring nonlocal indices
819: .seealso: ISGetTotalIndices(), ISGetNonlocalIndices(), ISRestoreTotalIndices()
820: @*/
821: PetscErrorCode  ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
822: {
826:   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.");
827:   return(0);
828: }

830: /*@
831:    ISGetNonlocalIS - Gather all nonlocal indices for this IS and present
832:                      them as another sequential index set.


835:    Collective on IS

837:    Input Parameter:
838: .  is - the index set

840:    Output Parameter:
841: .  complement - sequential IS with indices identical to the result of
842:                 ISGetNonlocalIndices()

844:    Level: intermediate

846:    Notes: complement represents the result of ISGetNonlocalIndices as an IS.
847:           Therefore scalability issues similar to ISGetNonlocalIndices apply.
848:           The resulting IS must be restored using ISRestoreNonlocalIS().

850:    Concepts: index sets^getting nonlocal indices
851: .seealso: ISGetNonlocalIndices(), ISRestoreNonlocalIndices(),  ISAllGather(), ISGetSize()
852: @*/
853: PetscErrorCode  ISGetNonlocalIS(IS is, IS *complement)
854: {

860:   /* Check if the complement exists already. */
861:   if (is->complement) {
862:     *complement = is->complement;
863:     PetscObjectReference((PetscObject)(is->complement));
864:   } else {
865:     PetscInt       N, n;
866:     const PetscInt *idx;
867:     ISGetSize(is, &N);
868:     ISGetLocalSize(is,&n);
869:     ISGetNonlocalIndices(is, &idx);
870:     ISCreateGeneral(PETSC_COMM_SELF, N-n,idx, PETSC_USE_POINTER, &(is->complement));
871:     PetscObjectReference((PetscObject)is->complement);
872:     *complement = is->complement;
873:   }
874:   return(0);
875: }


878: /*@
879:    ISRestoreNonlocalIS - Restore the IS obtained with ISGetNonlocalIS().

881:    Not collective.

883:    Input Parameter:
884: +  is         - the index set
885: -  complement - index set of is's nonlocal indices

887:    Level: intermediate


890:    Concepts: index sets^getting nonlocal indices
891:    Concepts: index sets^restoring nonlocal indices
892: .seealso: ISGetNonlocalIS(), ISGetNonlocalIndices(), ISRestoreNonlocalIndices()
893: @*/
894: PetscErrorCode  ISRestoreNonlocalIS(IS is, IS *complement)
895: {
897:   PetscInt       refcnt;

902:   if (*complement != is->complement) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Complement IS being restored was not obtained with ISGetNonlocalIS()");
903:   PetscObjectGetReference((PetscObject)(is->complement), &refcnt);
904:   if (refcnt <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Duplicate call to ISRestoreNonlocalIS() detected");
905:   PetscObjectDereference((PetscObject)(is->complement));
906:   return(0);
907: }

909: /*@C
910:    ISView - Displays an index set.

912:    Collective on IS

914:    Input Parameters:
915: +  is - the index set
916: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

918:    Level: intermediate

920: .seealso: PetscViewerASCIIOpen()
921: @*/
922: PetscErrorCode  ISView(IS is,PetscViewer viewer)
923: {

928:   if (!viewer) {
929:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is),&viewer);
930:   }

934:   PetscObjectPrintClassNamePrefixType((PetscObject)is,viewer);
935:   (*is->ops->view)(is,viewer);
936:   return(0);
937: }

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

942:   Collective on PetscViewer

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

948:   Level: intermediate

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

955:   Concepts: index set^loading from file

957: .seealso: PetscViewerBinaryOpen(), ISView(), MatLoad(), VecLoad()
958: @*/
959: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
960: {
961:   PetscBool      isbinary, ishdf5;

967:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
968:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
969:   if (!isbinary && !ishdf5) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
970:   if (!((PetscObject)is)->type_name) {ISSetType(is, ISGENERAL);}
971:   (*is->ops->load)(is, viewer);
972:   return(0);
973: }

975: /*@
976:    ISSort - Sorts the indices of an index set.

978:    Collective on IS

980:    Input Parameters:
981: .  is - the index set

983:    Level: intermediate

985:    Concepts: index sets^sorting
986:    Concepts: sorting^index set

988: .seealso: ISSortRemoveDups(), ISSorted()
989: @*/
990: PetscErrorCode  ISSort(IS is)
991: {

996:   (*is->ops->sort)(is);
997:   return(0);
998: }

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

1003:   Collective on IS

1005:   Input Parameters:
1006: . is - the index set

1008:   Level: intermediate

1010:   Concepts: index sets^sorting
1011:   Concepts: sorting^index set

1013: .seealso: ISSort(), ISSorted()
1014: @*/
1015: PetscErrorCode ISSortRemoveDups(IS is)
1016: {

1021:   (*is->ops->sortremovedups)(is);
1022:   return(0);
1023: }

1025: /*@
1026:    ISToGeneral - Converts an IS object of any type to ISGENERAL type

1028:    Collective on IS

1030:    Input Parameters:
1031: .  is - the index set

1033:    Level: intermediate

1035:    Concepts: index sets^sorting
1036:    Concepts: sorting^index set

1038: .seealso: ISSorted()
1039: @*/
1040: PetscErrorCode  ISToGeneral(IS is)
1041: {

1046:   if (is->ops->togeneral) {
1047:     (*is->ops->togeneral)(is);
1048:   } else SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Not written for this type %s",((PetscObject)is)->type_name);
1049:   return(0);
1050: }

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

1055:    Collective on IS

1057:    Input Parameter:
1058: .  is - the index set

1060:    Output Parameter:
1061: .  flg - output flag, either PETSC_TRUE if the index set is sorted,
1062:          or PETSC_FALSE otherwise.

1064:    Notes: For parallel IS objects this only indicates if the local part of the IS
1065:           is sorted. So some processors may return PETSC_TRUE while others may
1066:           return PETSC_FALSE.

1068:    Level: intermediate

1070: .seealso: ISSort(), ISSortRemoveDups()
1071: @*/
1072: PetscErrorCode  ISSorted(IS is,PetscBool  *flg)
1073: {

1079:   (*is->ops->sorted)(is,flg);
1080:   return(0);
1081: }

1083: /*@
1084:    ISDuplicate - Creates a duplicate copy of an index set.

1086:    Collective on IS

1088:    Input Parmeters:
1089: .  is - the index set

1091:    Output Parameters:
1092: .  isnew - the copy of the index set

1094:    Level: beginner

1096:    Concepts: index sets^duplicating

1098: .seealso: ISCreateGeneral(), ISCopy()
1099: @*/
1100: PetscErrorCode  ISDuplicate(IS is,IS *newIS)
1101: {

1107:   (*is->ops->duplicate)(is,newIS);
1108:   (*newIS)->isidentity = is->isidentity;
1109:   (*newIS)->isperm     = is->isperm;
1110:   return(0);
1111: }

1113: /*@
1114:    ISCopy - Copies an index set.

1116:    Collective on IS

1118:    Input Parmeters:
1119: .  is - the index set

1121:    Output Parameters:
1122: .  isy - the copy of the index set

1124:    Level: beginner

1126:    Concepts: index sets^copying

1128: .seealso: ISDuplicate()
1129: @*/
1130: PetscErrorCode  ISCopy(IS is,IS isy)
1131: {

1138:   if (is == isy) return(0);
1139:   (*is->ops->copy)(is,isy);
1140:   isy->isperm     = is->isperm;
1141:   isy->max        = is->max;
1142:   isy->min        = is->min;
1143:   isy->isidentity = is->isidentity;
1144:   return(0);
1145: }

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

1150:    Collective on IS and comm

1152:    Input Arguments:
1153: + is - index set
1154: . comm - communicator for new index set
1155: - mode - copy semantics, PETSC_USE_POINTER for no-copy if possible, otherwise PETSC_COPY_VALUES

1157:    Output Arguments:
1158: . newis - new IS on comm

1160:    Level: advanced

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

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

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

1170: .seealso: ISSplit()
1171: @*/
1172: PetscErrorCode  ISOnComm(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
1173: {
1175:   PetscMPIInt    match;

1180:   MPI_Comm_compare(PetscObjectComm((PetscObject)is),comm,&match);
1181:   if (mode != PETSC_COPY_VALUES && (match == MPI_IDENT || match == MPI_CONGRUENT)) {
1182:     PetscObjectReference((PetscObject)is);
1183:     *newis = is;
1184:   } else {
1185:     (*is->ops->oncomm)(is,comm,mode,newis);
1186:   }
1187:   return(0);
1188: }

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

1193:    Logicall Collective on IS

1195:    Input Arguments:
1196: + is - index set
1197: - bs - block size

1199:    Level: intermediate

1201: .seealso: ISGetBlockSize(), ISCreateBlock()
1202: @*/
1203: PetscErrorCode  ISSetBlockSize(IS is,PetscInt bs)
1204: {

1210:   if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_OUTOFRANGE,"Block size %D, must be positive",bs);
1211:   (*is->ops->setblocksize)(is,bs);
1212:   return(0);
1213: }

1215: /*@
1216:    ISGetBlockSize - Returns the number of elements in a block.

1218:    Not Collective

1220:    Input Parameter:
1221: .  is - the index set

1223:    Output Parameter:
1224: .  size - the number of elements in a block

1226:    Level: intermediate

1228:    Concepts: IS^block size
1229:    Concepts: index sets^block size

1231: .seealso: ISBlockGetSize(), ISGetSize(), ISCreateBlock(), ISSetBlockSize()
1232: @*/
1233: PetscErrorCode  ISGetBlockSize(IS is,PetscInt *size)
1234: {

1238:   PetscLayoutGetBlockSize(is->map, size);
1239:   return(0);
1240: }

1242: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[])
1243: {
1245:   PetscInt       len,i;
1246:   const PetscInt *ptr;

1249:   ISGetSize(is,&len);
1250:   ISGetIndices(is,&ptr);
1251:   for (i=0; i<len; i++) idx[i] = ptr[i];
1252:   ISRestoreIndices(is,&ptr);
1253:   return(0);
1254: }

1256: /*MC
1257:     ISGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1258:     The users should call ISRestoreIndicesF90() after having looked at the
1259:     indices.  The user should NOT change the indices.

1261:     Synopsis:
1262:     ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1264:     Not collective

1266:     Input Parameter:
1267: .   x - index set

1269:     Output Parameters:
1270: +   xx_v - the Fortran90 pointer to the array
1271: -   ierr - error code

1273:     Example of Usage:
1274: .vb
1275:     PetscScalar, pointer xx_v(:)
1276:     ....
1277:     call ISGetIndicesF90(x,xx_v,ierr)
1278:     a = xx_v(3)
1279:     call ISRestoreIndicesF90(x,xx_v,ierr)
1280: .ve

1282:     Notes:
1283:     Not yet supported for all F90 compilers.

1285:     Level: intermediate

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

1289:   Concepts: index sets^getting indices in f90
1290:   Concepts: indices of index set in f90

1292: M*/

1294: /*MC
1295:     ISRestoreIndicesF90 - Restores an index set to a usable state after
1296:     a call to ISGetIndicesF90().

1298:     Synopsis:
1299:     ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1301:     Not collective

1303:     Input Parameters:
1304: .   x - index set
1305: .   xx_v - the Fortran90 pointer to the array

1307:     Output Parameter:
1308: .   ierr - error code


1311:     Example of Usage:
1312: .vb
1313:     PetscScalar, pointer xx_v(:)
1314:     ....
1315:     call ISGetIndicesF90(x,xx_v,ierr)
1316:     a = xx_v(3)
1317:     call ISRestoreIndicesF90(x,xx_v,ierr)
1318: .ve

1320:     Notes:
1321:     Not yet supported for all F90 compilers.

1323:     Level: intermediate

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

1327: M*/

1329: /*MC
1330:     ISBlockGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1331:     The users should call ISBlockRestoreIndicesF90() after having looked at the
1332:     indices.  The user should NOT change the indices.

1334:     Synopsis:
1335:     ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1337:     Not collective

1339:     Input Parameter:
1340: .   x - index set

1342:     Output Parameters:
1343: +   xx_v - the Fortran90 pointer to the array
1344: -   ierr - error code
1345:     Example of Usage:
1346: .vb
1347:     PetscScalar, pointer xx_v(:)
1348:     ....
1349:     call ISBlockGetIndicesF90(x,xx_v,ierr)
1350:     a = xx_v(3)
1351:     call ISBlockRestoreIndicesF90(x,xx_v,ierr)
1352: .ve

1354:     Notes:
1355:     Not yet supported for all F90 compilers

1357:     Level: intermediate

1359: .seealso:  ISBlockRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices(),
1360:            ISRestoreIndices()

1362:   Concepts: index sets^getting block indices in f90
1363:   Concepts: indices of index set in f90
1364:   Concepts: block^ indices of index set in f90

1366: M*/

1368: /*MC
1369:     ISBlockRestoreIndicesF90 - Restores an index set to a usable state after
1370:     a call to ISBlockGetIndicesF90().

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

1375:     Not Collective

1377:     Input Parameters:
1378: +   x - index set
1379: -   xx_v - the Fortran90 pointer to the array

1381:     Output Parameter:
1382: .   ierr - error code

1384:     Example of Usage:
1385: .vb
1386:     PetscScalar, pointer xx_v(:)
1387:     ....
1388:     call ISBlockGetIndicesF90(x,xx_v,ierr)
1389:     a = xx_v(3)
1390:     call ISBlockRestoreIndicesF90(x,xx_v,ierr)
1391: .ve

1393:     Notes:
1394:     Not yet supported for all F90 compilers

1396:     Level: intermediate

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

1400: M*/