Actual source code: index.c

petsc-3.12.1 2019-10-22
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_mult,&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;

 64:   PetscLayoutCreate(PetscObjectComm((PetscObject)subset),&map);
 65:   PetscLayoutSetBlockSize(map,1);
 66:   PetscLayoutSetSize(map,N_n);
 67:   PetscLayoutSetUp(map);
 68:   PetscLayoutGetLocalSize(map,&Nl);

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

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

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

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

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

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

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

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

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

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

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

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

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


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

177:    Collective on IS

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

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

186:    Level: intermediate

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

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


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

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


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

252:    Collective on IS

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

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

260:    Level: intermediate


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

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

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

283:    Logically Collective on IS

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

288:    Level: intermediate


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

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

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

307:    Not Collective

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

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

318:    Level: developer

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

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

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

343:    Logically Collective on IS

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

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

351:    Level: intermediate


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

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

368:    Logically Collective on IS

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

373:    Level: intermediate


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

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

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

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

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

415:    Collective on IS

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

420:    Level: beginner

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

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

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

453:    Collective on IS

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

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

463:    Level: intermediate

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

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

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

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

490:    Not Collective

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

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

498:    Level: beginner


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

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

514:    Not Collective

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

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

522:    Level: beginner

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

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

539:    Not Collective

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

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

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



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

572:    Level: intermediate


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

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

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

591:    Not Collective

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

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

600:    Level: intermediate

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


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

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

621:   Not Collective

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

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

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

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

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

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

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

669:    Not Collective

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

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

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

692:    Level: intermediate

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

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

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

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


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

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

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

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

747:    Collective on IS

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

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

756:    Level: intermediate

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

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

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

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

791:    Not Collective.

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

797:    Level: intermediate

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

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

821:    Collective on IS

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

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

832:    Level: intermediate

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

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

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

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

869:    Not Collective.

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

875:    Level: intermediate

877: .seealso: ISGetTotalIndices(), ISGetNonlocalIndices(), ISRestoreTotalIndices()
878: @*/
879: PetscErrorCode  ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
880: {
884:   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.");
885:   return(0);
886: }

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


893:    Collective on IS

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

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

902:    Level: intermediate

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

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

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


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

939:    Not collective.

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

945:    Level: intermediate


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

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

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

968:    Collective on IS

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

974:    Level: intermediate

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

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

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

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

998:   Collective on PetscViewer

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

1004:   Level: intermediate

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

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

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

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

1032:    Collective on IS

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

1037:    Level: intermediate


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

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

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

1055:   Collective on IS

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

1060:   Level: intermediate


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

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

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

1078:    Collective on IS

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

1083:    Level: intermediate


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

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

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

1103:    Collective on IS

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

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

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

1117:    Level: intermediate

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

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

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

1135:    Collective on IS

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

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

1143:    Level: beginner

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

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

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

1163:    Collective on IS

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

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

1171:    Level: beginner

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

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

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

1195:    Collective on IS

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

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

1205:    Level: advanced

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

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

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

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

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

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

1238:    Logicall Collective on IS

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

1244:    Level: intermediate

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

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

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

1263:    Not Collective

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

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

1271:    Level: intermediate


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

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

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

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

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

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

1307:     Not collective

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

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

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

1325:     Level: intermediate

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


1330: M*/

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

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

1339:     Not collective

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

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


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

1358:     Level: intermediate

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

1362: M*/

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

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

1372:     Not collective

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

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

1389:     Level: intermediate

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

1394: M*/

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

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

1403:     Not Collective

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

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

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

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

1424:     Level: intermediate

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

1428: M*/