Actual source code: index.c

petsc-master 2016-07-26
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>      /*I "petscis.h" I*/
  6: #include <petscviewer.h>
  7: #include <petscsf.h>

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

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

 17:    Collective on IS

 19:    Input Parmeters:
 20: .  subset - the index set
 21: .  subset_mult - the multiplcity of each entry in subset (optional, can be NULL)

 23:    Output Parameters:
 24: .  N - the maximum entry of the new IS
 25: .  subset_n - the new IS

 27:    Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

165: /*@
166:    ISIdentity - Determines whether index set is the identity mapping.

168:    Collective on IS

170:    Input Parmeters:
171: .  is - the index set

173:    Output Parameters:
174: .  ident - PETSC_TRUE if an identity, else PETSC_FALSE

176:    Level: intermediate

178:    Concepts: identity mapping
179:    Concepts: index sets^is identity

181: .seealso: ISSetIdentity()
182: @*/
183: PetscErrorCode  ISIdentity(IS is,PetscBool  *ident)
184: {

190:   *ident = is->isidentity;
191:   if (*ident) return(0);
192:   if (is->ops->identity) {
193:     (*is->ops->identity)(is,ident);
194:   }
195:   return(0);
196: }

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

203:    Logically Collective on IS

205:    Input Parmeters:
206: .  is - the index set

208:    Level: intermediate

210:    Concepts: identity mapping
211:    Concepts: index sets^is identity

213: .seealso: ISIdentity()
214: @*/
215: PetscErrorCode  ISSetIdentity(IS is)
216: {

221:   is->isidentity = PETSC_TRUE;
222:   ISSetPermutation(is);
223:   return(0);
224: }

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

231:    Not Collective

233:    Input Parmeters:
234: +  is - the index set
235: .  gstart - global start
236: .  gend - global end

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

242:    Level: developer

244:    Concepts: index sets^is contiguous

246: .seealso: ISGetLocalSize(), VecGetOwnershipRange()
247: @*/
248: PetscErrorCode  ISContiguousLocal(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
249: {

256:   if (is->ops->contiguous) {
257:     (*is->ops->contiguous)(is,gstart,gend,start,contig);
258:   } else {
259:     *start  = -1;
260:     *contig = PETSC_FALSE;
261:   }
262:   return(0);
263: }

267: /*@
268:    ISPermutation - PETSC_TRUE or PETSC_FALSE depending on whether the
269:    index set has been declared to be a permutation.

271:    Logically Collective on IS

273:    Input Parmeters:
274: .  is - the index set

276:    Output Parameters:
277: .  perm - PETSC_TRUE if a permutation, else PETSC_FALSE

279:    Level: intermediate

281:   Concepts: permutation
282:   Concepts: index sets^is permutation

284: .seealso: ISSetPermutation()
285: @*/
286: PetscErrorCode  ISPermutation(IS is,PetscBool  *perm)
287: {
291:   *perm = (PetscBool) is->isperm;
292:   return(0);
293: }

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

300:    Logically Collective on IS

302:    Input Parmeters:
303: .  is - the index set

305:    Level: intermediate

307:   Concepts: permutation
308:   Concepts: index sets^permutation

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

313: .seealso: ISPermutation()
314: @*/
315: PetscErrorCode  ISSetPermutation(IS is)
316: {
319: #if defined(PETSC_USE_DEBUG)
320:   {
321:     PetscMPIInt    size;

324:     MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
325:     if (size == 1) {
326:       PetscInt       i,n,*idx;
327:       const PetscInt *iidx;

329:       ISGetSize(is,&n);
330:       PetscMalloc1(n,&idx);
331:       ISGetIndices(is,&iidx);
332:       PetscMemcpy(idx,iidx,n*sizeof(PetscInt));
333:       PetscSortInt(n,idx);
334:       for (i=0; i<n; i++) {
335:         if (idx[i] != i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Index set is not a permutation");
336:       }
337:       PetscFree(idx);
338:       ISRestoreIndices(is,&iidx);
339:     }
340:   }
341: #endif
342:   is->isperm = PETSC_TRUE;
343:   return(0);
344: }

348: /*@
349:    ISDestroy - Destroys an index set.

351:    Collective on IS

353:    Input Parameters:
354: .  is - the index set

356:    Level: beginner

358: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlocked()
359: @*/
360: PetscErrorCode  ISDestroy(IS *is)
361: {

365:   if (!*is) return(0);
367:   if (--((PetscObject)(*is))->refct > 0) {*is = 0; return(0);}
368:   if ((*is)->complement) {
369:     PetscInt refcnt;
370:     PetscObjectGetReference((PetscObject)((*is)->complement), &refcnt);
371:     if (refcnt > 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Nonlocal IS has not been restored");
372:     ISDestroy(&(*is)->complement);
373:   }
374:   if ((*is)->ops->destroy) {
375:     (*(*is)->ops->destroy)(*is);
376:   }
377:   PetscLayoutDestroy(&(*is)->map);
378:   /* Destroy local representations of offproc data. */
379:   PetscFree((*is)->total);
380:   PetscFree((*is)->nonlocal);
381:   PetscHeaderDestroy(is);
382:   return(0);
383: }

387: /*@
388:    ISInvertPermutation - Creates a new permutation that is the inverse of
389:                          a given permutation.

391:    Collective on IS

393:    Input Parameter:
394: +  is - the index set
395: -  nlocal - number of indices on this processor in result (ignored for 1 proccessor) or
396:             use PETSC_DECIDE

398:    Output Parameter:
399: .  isout - the inverse permutation

401:    Level: intermediate

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

406:    Concepts: inverse permutation
407:    Concepts: permutation^inverse
408:    Concepts: index sets^inverting
409: @*/
410: PetscErrorCode  ISInvertPermutation(IS is,PetscInt nlocal,IS *isout)
411: {

417:   if (!is->isperm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a permutation, must call ISSetPermutation() on the IS first");
418:   if (is->isidentity) {
419:     ISDuplicate(is,isout);
420:   } else {
421:     (*is->ops->invertpermutation)(is,nlocal,isout);
422:     ISSetPermutation(*isout);
423:   }
424:   return(0);
425: }

429: /*@
430:    ISGetSize - Returns the global length of an index set.

432:    Not Collective

434:    Input Parameter:
435: .  is - the index set

437:    Output Parameter:
438: .  size - the global size

440:    Level: beginner

442:    Concepts: size^of index set
443:    Concepts: index sets^size

445: @*/
446: PetscErrorCode  ISGetSize(IS is,PetscInt *size)
447: {

453:   (*is->ops->getsize)(is,size);
454:   return(0);
455: }

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

462:    Not Collective

464:    Input Parameter:
465: .  is - the index set

467:    Output Parameter:
468: .  size - the local size

470:    Level: beginner

472:    Concepts: size^of index set
473:    Concepts: local size^of index set
474:    Concepts: index sets^local size

476: @*/
477: PetscErrorCode  ISGetLocalSize(IS is,PetscInt *size)
478: {

484:   (*is->ops->getlocalsize)(is,size);
485:   return(0);
486: }

490: /*@C
491:    ISGetIndices - Returns a pointer to the indices.  The user should call
492:    ISRestoreIndices() after having looked at the indices.  The user should
493:    NOT change the indices.

495:    Not Collective

497:    Input Parameter:
498: .  is - the index set

500:    Output Parameter:
501: .  ptr - the location to put the pointer to the indices

503:    Fortran Note:
504:    This routine is used differently from Fortran
505: $    IS          is
506: $    integer     is_array(1)
507: $    PetscOffset i_is
508: $    int         ierr
509: $       call ISGetIndices(is,is_array,i_is,ierr)
510: $
511: $   Access first local entry in list
512: $      value = is_array(i_is + 1)
513: $
514: $      ...... other code
515: $       call ISRestoreIndices(is,is_array,i_is,ierr)

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

520:    Level: intermediate

522:    Concepts: index sets^getting indices
523:    Concepts: indices of index set

525: .seealso: ISRestoreIndices(), ISGetIndicesF90()
526: @*/
527: PetscErrorCode  ISGetIndices(IS is,const PetscInt *ptr[])
528: {

534:   (*is->ops->getindices)(is,ptr);
535:   return(0);
536: }

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

543:    Not Collective

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

548:    Output Parameter:
549: +   min - the minimum value
550: -   max - the maximum value

552:    Level: intermediate

554:    Concepts: index sets^getting indices
555:    Concepts: indices of index set

557: .seealso: ISGetIndices(), ISRestoreIndices(), ISGetIndicesF90()
558: @*/
559: PetscErrorCode  ISGetMinMax(IS is,PetscInt *min,PetscInt *max)
560: {
563:   if (min) *min = is->min;
564:   if (max) *max = is->max;
565:   return(0);
566: }

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

573:   Not Collective

575:   Input Parameter:
576: + is - the index set
577: - key - the search key

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

582:   Level: intermediate
583: @*/
584: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
585: {

589:   if (is->ops->locate) {
590:     (*is->ops->locate)(is,key,location);
591:   } else {
592:     PetscInt       numIdx;
593:     PetscBool      sorted;
594:     const PetscInt *idx;

596:     ISGetLocalSize(is,&numIdx);
597:     ISGetIndices(is,&idx);
598:     ISSorted(is,&sorted);
599:     if (sorted) {
600:       PetscFindInt(key,numIdx,idx,location);
601:     } else {
602:       PetscInt i;

604:       *location = -1;
605:       for (i = 0; i < numIdx; i++) {
606:         if (idx[i] == key) {
607:           *location = i;
608:           break;
609:         }
610:       }
611:     }
612:     ISRestoreIndices(is,&idx);
613:   }
614:   return(0);
615: }

619: /*@C
620:    ISRestoreIndices - Restores an index set to a usable state after a call
621:                       to ISGetIndices().

623:    Not Collective

625:    Input Parameters:
626: +  is - the index set
627: -  ptr - the pointer obtained by ISGetIndices()

629:    Fortran Note:
630:    This routine is used differently from Fortran
631: $    IS          is
632: $    integer     is_array(1)
633: $    PetscOffset i_is
634: $    int         ierr
635: $       call ISGetIndices(is,is_array,i_is,ierr)
636: $
637: $   Access first local entry in list
638: $      value = is_array(i_is + 1)
639: $
640: $      ...... other code
641: $       call ISRestoreIndices(is,is_array,i_is,ierr)

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

646:    Level: intermediate

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

651: .seealso: ISGetIndices(), ISRestoreIndicesF90()
652: @*/
653: PetscErrorCode  ISRestoreIndices(IS is,const PetscInt *ptr[])
654: {

660:   if (is->ops->restoreindices) {
661:     (*is->ops->restoreindices)(is,ptr);
662:   }
663:   return(0);
664: }

668: static PetscErrorCode ISGatherTotal_Private(IS is)
669: {
671:   PetscInt       i,n,N;
672:   const PetscInt *lindices;
673:   MPI_Comm       comm;
674:   PetscMPIInt    rank,size,*sizes = NULL,*offsets = NULL,nn;


679:   PetscObjectGetComm((PetscObject)is,&comm);
680:   MPI_Comm_size(comm,&size);
681:   MPI_Comm_rank(comm,&rank);
682:   ISGetLocalSize(is,&n);
683:   PetscMalloc2(size,&sizes,size,&offsets);

685:   PetscMPIIntCast(n,&nn);
686:   MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
687:   offsets[0] = 0;
688:   for (i=1; i<size; ++i) offsets[i] = offsets[i-1] + sizes[i-1];
689:   N = offsets[size-1] + sizes[size-1];

691:   PetscMalloc1(N,&(is->total));
692:   ISGetIndices(is,&lindices);
693:   MPI_Allgatherv((void*)lindices,nn,MPIU_INT,is->total,sizes,offsets,MPIU_INT,comm);
694:   ISRestoreIndices(is,&lindices);
695:   is->local_offset = offsets[rank];
696:   PetscFree2(sizes,offsets);
697:   return(0);
698: }

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

705:    Collective on IS

707:    Input Parameter:
708: .  is - the index set

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

714:    Level: intermediate

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

723:    Concepts: index sets^getting nonlocal indices
724: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices(), ISGetSize()
725: @*/
726: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
727: {
729:   PetscMPIInt    size;

734:   MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
735:   if (size == 1) {
736:     (*is->ops->getindices)(is,indices);
737:   } else {
738:     if (!is->total) {
739:       ISGatherTotal_Private(is);
740:     }
741:     *indices = is->total;
742:   }
743:   return(0);
744: }

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

751:    Not Collective.

753:    Input Parameter:
754: +  is - the index set
755: -  indices - index array; must be the array obtained with ISGetTotalIndices()

757:    Level: intermediate

759:    Concepts: index sets^getting nonlocal indices
760:    Concepts: index sets^restoring nonlocal indices
761: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices()
762: @*/
763: PetscErrorCode  ISRestoreTotalIndices(IS is, const PetscInt *indices[])
764: {
766:   PetscMPIInt    size;

771:   MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
772:   if (size == 1) {
773:     (*is->ops->restoreindices)(is,indices);
774:   } else {
775:     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.");
776:   }
777:   return(0);
778: }
781: /*@C
782:    ISGetNonlocalIndices - Retrieve an array of indices from remote processors
783:                        in this communicator.

785:    Collective on IS

787:    Input Parameter:
788: .  is - the index set

790:    Output Parameter:
791: .  indices - indices with rank 0 indices first, and so on,  omitting
792:              the current rank.  Total number of indices is the difference
793:              total and local, obtained with ISGetSize() and ISGetLocalSize(),
794:              respectively.

796:    Level: intermediate

798:    Notes: restore the indices using ISRestoreNonlocalIndices().
799:           The same scalability considerations as those for ISGetTotalIndices
800:           apply here.

802:    Concepts: index sets^getting nonlocal indices
803: .seealso: ISGetTotalIndices(), ISRestoreNonlocalIndices(), ISGetSize(), ISGetLocalSize().
804: @*/
805: PetscErrorCode  ISGetNonlocalIndices(IS is, const PetscInt *indices[])
806: {
808:   PetscMPIInt    size;
809:   PetscInt       n, N;

814:   MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
815:   if (size == 1) *indices = NULL;
816:   else {
817:     if (!is->total) {
818:       ISGatherTotal_Private(is);
819:     }
820:     ISGetLocalSize(is,&n);
821:     ISGetSize(is,&N);
822:     PetscMalloc1(N-n, &(is->nonlocal));
823:     PetscMemcpy(is->nonlocal, is->total, sizeof(PetscInt)*is->local_offset);
824:     PetscMemcpy(is->nonlocal+is->local_offset, is->total+is->local_offset+n, sizeof(PetscInt)*(N - is->local_offset - n));
825:     *indices = is->nonlocal;
826:   }
827:   return(0);
828: }

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

835:    Not Collective.

837:    Input Parameter:
838: +  is - the index set
839: -  indices - index array; must be the array obtained with ISGetNonlocalIndices()

841:    Level: intermediate

843:    Concepts: index sets^getting nonlocal indices
844:    Concepts: index sets^restoring nonlocal indices
845: .seealso: ISGetTotalIndices(), ISGetNonlocalIndices(), ISRestoreTotalIndices()
846: @*/
847: PetscErrorCode  ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
848: {
852:   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.");
853:   return(0);
854: }

858: /*@
859:    ISGetNonlocalIS - Gather all nonlocal indices for this IS and present
860:                      them as another sequential index set.


863:    Collective on IS

865:    Input Parameter:
866: .  is - the index set

868:    Output Parameter:
869: .  complement - sequential IS with indices identical to the result of
870:                 ISGetNonlocalIndices()

872:    Level: intermediate

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

878:    Concepts: index sets^getting nonlocal indices
879: .seealso: ISGetNonlocalIndices(), ISRestoreNonlocalIndices(),  ISAllGather(), ISGetSize()
880: @*/
881: PetscErrorCode  ISGetNonlocalIS(IS is, IS *complement)
882: {

888:   /* Check if the complement exists already. */
889:   if (is->complement) {
890:     *complement = is->complement;
891:     PetscObjectReference((PetscObject)(is->complement));
892:   } else {
893:     PetscInt       N, n;
894:     const PetscInt *idx;
895:     ISGetSize(is, &N);
896:     ISGetLocalSize(is,&n);
897:     ISGetNonlocalIndices(is, &idx);
898:     ISCreateGeneral(PETSC_COMM_SELF, N-n,idx, PETSC_USE_POINTER, &(is->complement));
899:     PetscObjectReference((PetscObject)is->complement);
900:     *complement = is->complement;
901:   }
902:   return(0);
903: }


908: /*@
909:    ISRestoreNonlocalIS - Restore the IS obtained with ISGetNonlocalIS().

911:    Not collective.

913:    Input Parameter:
914: +  is         - the index set
915: -  complement - index set of is's nonlocal indices

917:    Level: intermediate


920:    Concepts: index sets^getting nonlocal indices
921:    Concepts: index sets^restoring nonlocal indices
922: .seealso: ISGetNonlocalIS(), ISGetNonlocalIndices(), ISRestoreNonlocalIndices()
923: @*/
924: PetscErrorCode  ISRestoreNonlocalIS(IS is, IS *complement)
925: {
927:   PetscInt       refcnt;

932:   if (*complement != is->complement) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Complement IS being restored was not obtained with ISGetNonlocalIS()");
933:   PetscObjectGetReference((PetscObject)(is->complement), &refcnt);
934:   if (refcnt <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Duplicate call to ISRestoreNonlocalIS() detected");
935:   PetscObjectDereference((PetscObject)(is->complement));
936:   return(0);
937: }

941: /*@C
942:    ISView - Displays an index set.

944:    Collective on IS

946:    Input Parameters:
947: +  is - the index set
948: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

950:    Level: intermediate

952: .seealso: PetscViewerASCIIOpen()
953: @*/
954: PetscErrorCode  ISView(IS is,PetscViewer viewer)
955: {

960:   if (!viewer) {
961:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is),&viewer);
962:   }

966:   PetscObjectPrintClassNamePrefixType((PetscObject)is,viewer);
967:   (*is->ops->view)(is,viewer);
968:   return(0);
969: }

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

976:   Collective on PetscViewer

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

982:   Level: intermediate

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

989:   Concepts: index set^loading from file

991: .seealso: PetscViewerBinaryOpen(), ISView(), MatLoad(), VecLoad()
992: @*/
993: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
994: {
995:   PetscBool      isbinary, ishdf5;

1001:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
1002:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
1003:   if (!isbinary && !ishdf5) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1004:   if (!((PetscObject)is)->type_name) {ISSetType(is, ISGENERAL);}
1005:   (*is->ops->load)(is, viewer);
1006:   return(0);
1007: }

1011: /*@
1012:    ISSort - Sorts the indices of an index set.

1014:    Collective on IS

1016:    Input Parameters:
1017: .  is - the index set

1019:    Level: intermediate

1021:    Concepts: index sets^sorting
1022:    Concepts: sorting^index set

1024: .seealso: ISSortRemoveDups(), ISSorted()
1025: @*/
1026: PetscErrorCode  ISSort(IS is)
1027: {

1032:   (*is->ops->sort)(is);
1033:   return(0);
1034: }

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

1041:   Collective on IS

1043:   Input Parameters:
1044: . is - the index set

1046:   Level: intermediate

1048:   Concepts: index sets^sorting
1049:   Concepts: sorting^index set

1051: .seealso: ISSort(), ISSorted()
1052: @*/
1053: PetscErrorCode ISSortRemoveDups(IS is)
1054: {

1059:   (*is->ops->sortremovedups)(is);
1060:   return(0);
1061: }

1065: /*@
1066:    ISToGeneral - Converts an IS object of any type to ISGENERAL type

1068:    Collective on IS

1070:    Input Parameters:
1071: .  is - the index set

1073:    Level: intermediate

1075:    Concepts: index sets^sorting
1076:    Concepts: sorting^index set

1078: .seealso: ISSorted()
1079: @*/
1080: PetscErrorCode  ISToGeneral(IS is)
1081: {

1086:   if (is->ops->togeneral) {
1087:     (*is->ops->togeneral)(is);
1088:   } else SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Not written for this type %s",((PetscObject)is)->type_name);
1089:   return(0);
1090: }

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

1097:    Collective on IS

1099:    Input Parameter:
1100: .  is - the index set

1102:    Output Parameter:
1103: .  flg - output flag, either PETSC_TRUE if the index set is sorted,
1104:          or PETSC_FALSE otherwise.

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

1110:    Level: intermediate

1112: .seealso: ISSort(), ISSortRemoveDups()
1113: @*/
1114: PetscErrorCode  ISSorted(IS is,PetscBool  *flg)
1115: {

1121:   (*is->ops->sorted)(is,flg);
1122:   return(0);
1123: }

1127: /*@
1128:    ISDuplicate - Creates a duplicate copy of an index set.

1130:    Collective on IS

1132:    Input Parmeters:
1133: .  is - the index set

1135:    Output Parameters:
1136: .  isnew - the copy of the index set

1138:    Level: beginner

1140:    Concepts: index sets^duplicating

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

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

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

1162:    Collective on IS

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

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

1170:    Level: beginner

1172:    Concepts: index sets^copying

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

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

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 and comm

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: }

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

1243:    Logicall Collective on IS

1245:    Input Arguments:
1246: + is - index set
1247: - bs - block size

1249:    Level: intermediate

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

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

1267: /*@
1268:    ISGetBlockSize - Returns the number of elements in a block.

1270:    Not Collective

1272:    Input Parameter:
1273: .  is - the index set

1275:    Output Parameter:
1276: .  size - the number of elements in a block

1278:    Level: intermediate

1280:    Concepts: IS^block size
1281:    Concepts: index sets^block size

1283: .seealso: ISBlockGetSize(), ISGetSize(), ISCreateBlock(), ISSetBlockSize()
1284: @*/
1285: PetscErrorCode  ISGetBlockSize(IS is,PetscInt *size)
1286: {

1290:   PetscLayoutGetBlockSize(is->map, size);
1291:   return(0);
1292: }

1296: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[])
1297: {
1299:   PetscInt       len,i;
1300:   const PetscInt *ptr;

1303:   ISGetSize(is,&len);
1304:   ISGetIndices(is,&ptr);
1305:   for (i=0; i<len; i++) idx[i] = ptr[i];
1306:   ISRestoreIndices(is,&ptr);
1307:   return(0);
1308: }

1310: /*MC
1311:     ISGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1312:     The users should call ISRestoreIndicesF90() after having looked at the
1313:     indices.  The user should NOT change the indices.

1315:     Synopsis:
1316:     ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1318:     Not collective

1320:     Input Parameter:
1321: .   x - index set

1323:     Output Parameters:
1324: +   xx_v - the Fortran90 pointer to the array
1325: -   ierr - error code

1327:     Example of Usage:
1328: .vb
1329:     PetscScalar, pointer xx_v(:)
1330:     ....
1331:     call ISGetIndicesF90(x,xx_v,ierr)
1332:     a = xx_v(3)
1333:     call ISRestoreIndicesF90(x,xx_v,ierr)
1334: .ve

1336:     Notes:
1337:     Not yet supported for all F90 compilers.

1339:     Level: intermediate

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

1343:   Concepts: index sets^getting indices in f90
1344:   Concepts: indices of index set in f90

1346: M*/

1348: /*MC
1349:     ISRestoreIndicesF90 - Restores an index set to a usable state after
1350:     a call to ISGetIndicesF90().

1352:     Synopsis:
1353:     ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1355:     Not collective

1357:     Input Parameters:
1358: .   x - index set
1359: .   xx_v - the Fortran90 pointer to the array

1361:     Output Parameter:
1362: .   ierr - error code


1365:     Example of Usage:
1366: .vb
1367:     PetscScalar, pointer xx_v(:)
1368:     ....
1369:     call ISGetIndicesF90(x,xx_v,ierr)
1370:     a = xx_v(3)
1371:     call ISRestoreIndicesF90(x,xx_v,ierr)
1372: .ve

1374:     Notes:
1375:     Not yet supported for all F90 compilers.

1377:     Level: intermediate

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

1381: M*/

1383: /*MC
1384:     ISBlockGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1385:     The users should call ISBlockRestoreIndicesF90() after having looked at the
1386:     indices.  The user should NOT change the indices.

1388:     Synopsis:
1389:     ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1391:     Not collective

1393:     Input Parameter:
1394: .   x - index set

1396:     Output Parameters:
1397: +   xx_v - the Fortran90 pointer to the array
1398: -   ierr - error code
1399:     Example of Usage:
1400: .vb
1401:     PetscScalar, pointer xx_v(:)
1402:     ....
1403:     call ISBlockGetIndicesF90(x,xx_v,ierr)
1404:     a = xx_v(3)
1405:     call ISBlockRestoreIndicesF90(x,xx_v,ierr)
1406: .ve

1408:     Notes:
1409:     Not yet supported for all F90 compilers

1411:     Level: intermediate

1413: .seealso:  ISBlockRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices(),
1414:            ISRestoreIndices()

1416:   Concepts: index sets^getting block indices in f90
1417:   Concepts: indices of index set in f90
1418:   Concepts: block^ indices of index set in f90

1420: M*/

1422: /*MC
1423:     ISBlockRestoreIndicesF90 - Restores an index set to a usable state after
1424:     a call to ISBlockGetIndicesF90().

1426:     Synopsis:
1427:     ISBlockRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1429:     Not Collective

1431:     Input Parameters:
1432: +   x - index set
1433: -   xx_v - the Fortran90 pointer to the array

1435:     Output Parameter:
1436: .   ierr - error code

1438:     Example of Usage:
1439: .vb
1440:     PetscScalar, pointer xx_v(:)
1441:     ....
1442:     call ISBlockGetIndicesF90(x,xx_v,ierr)
1443:     a = xx_v(3)
1444:     call ISBlockRestoreIndicesF90(x,xx_v,ierr)
1445: .ve

1447:     Notes:
1448:     Not yet supported for all F90 compilers

1450:     Level: intermediate

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

1454: M*/