Actual source code: index.c

petsc-master 2019-05-18
Report Typos and Errors

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

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

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

 15:    Collective on IS

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

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

 25:    Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

176:    Collective on IS

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

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

185:    Level: intermediate

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

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


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

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


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

251:    Collective on IS

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

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

259:    Level: intermediate

261:    Concepts: identity mapping
262:    Concepts: index sets^is identity

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

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

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

284:    Logically Collective on IS

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

289:    Level: intermediate

291:    Concepts: identity mapping
292:    Concepts: index sets^is identity

294: .seealso: ISIdentity()
295: @*/
296: PetscErrorCode  ISSetIdentity(IS is)
297: {

302:   is->isidentity = PETSC_TRUE;
303:   ISSetPermutation(is);
304:   return(0);
305: }

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

310:    Not Collective

312:    Input Parmeters:
313: +  is - the index set
314: .  gstart - global start
315: .  gend - global end

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

321:    Level: developer

323:    Concepts: index sets^is contiguous

325: .seealso: ISGetLocalSize(), VecGetOwnershipRange()
326: @*/
327: PetscErrorCode  ISContiguousLocal(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
328: {

335:   if (is->ops->contiguous) {
336:     (*is->ops->contiguous)(is,gstart,gend,start,contig);
337:   } else {
338:     *start  = -1;
339:     *contig = PETSC_FALSE;
340:   }
341:   return(0);
342: }

344: /*@
345:    ISPermutation - PETSC_TRUE or PETSC_FALSE depending on whether the
346:    index set has been declared to be a permutation.

348:    Logically Collective on IS

350:    Input Parmeters:
351: .  is - the index set

353:    Output Parameters:
354: .  perm - PETSC_TRUE if a permutation, else PETSC_FALSE

356:    Level: intermediate

358:   Concepts: permutation
359:   Concepts: index sets^is permutation

361: .seealso: ISSetPermutation()
362: @*/
363: PetscErrorCode  ISPermutation(IS is,PetscBool  *perm)
364: {
368:   *perm = (PetscBool) is->isperm;
369:   return(0);
370: }

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

375:    Logically Collective on IS

377:    Input Parmeters:
378: .  is - the index set

380:    Level: intermediate

382:   Concepts: permutation
383:   Concepts: index sets^permutation

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

388: .seealso: ISPermutation()
389: @*/
390: PetscErrorCode  ISSetPermutation(IS is)
391: {
394: #if defined(PETSC_USE_DEBUG)
395:   {
396:     PetscMPIInt    size;

399:     MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
400:     if (size == 1) {
401:       PetscInt       i,n,*idx;
402:       const PetscInt *iidx;

404:       ISGetSize(is,&n);
405:       PetscMalloc1(n,&idx);
406:       ISGetIndices(is,&iidx);
407:       PetscMemcpy(idx,iidx,n*sizeof(PetscInt));
408:       PetscSortInt(n,idx);
409:       for (i=0; i<n; i++) {
410:         if (idx[i] != i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Index set is not a permutation");
411:       }
412:       PetscFree(idx);
413:       ISRestoreIndices(is,&iidx);
414:     }
415:   }
416: #endif
417:   is->isperm = PETSC_TRUE;
418:   return(0);
419: }

421: /*@
422:    ISDestroy - Destroys an index set.

424:    Collective on IS

426:    Input Parameters:
427: .  is - the index set

429:    Level: beginner

431: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlocked()
432: @*/
433: PetscErrorCode  ISDestroy(IS *is)
434: {

438:   if (!*is) return(0);
440:   if (--((PetscObject)(*is))->refct > 0) {*is = 0; return(0);}
441:   if ((*is)->complement) {
442:     PetscInt refcnt;
443:     PetscObjectGetReference((PetscObject)((*is)->complement), &refcnt);
444:     if (refcnt > 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Nonlocal IS has not been restored");
445:     ISDestroy(&(*is)->complement);
446:   }
447:   if ((*is)->ops->destroy) {
448:     (*(*is)->ops->destroy)(*is);
449:   }
450:   PetscLayoutDestroy(&(*is)->map);
451:   /* Destroy local representations of offproc data. */
452:   PetscFree((*is)->total);
453:   PetscFree((*is)->nonlocal);
454:   PetscHeaderDestroy(is);
455:   return(0);
456: }

458: /*@
459:    ISInvertPermutation - Creates a new permutation that is the inverse of
460:                          a given permutation.

462:    Collective on IS

464:    Input Parameter:
465: +  is - the index set
466: -  nlocal - number of indices on this processor in result (ignored for 1 proccessor) or
467:             use PETSC_DECIDE

469:    Output Parameter:
470: .  isout - the inverse permutation

472:    Level: intermediate

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

478:    Concepts: inverse permutation
479:    Concepts: permutation^inverse
480:    Concepts: index sets^inverting
481: @*/
482: PetscErrorCode  ISInvertPermutation(IS is,PetscInt nlocal,IS *isout)
483: {

489:   if (!is->isperm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a permutation, must call ISSetPermutation() on the IS first");
490:   if (is->isidentity) {
491:     ISDuplicate(is,isout);
492:   } else {
493:     (*is->ops->invertpermutation)(is,nlocal,isout);
494:     ISSetPermutation(*isout);
495:   }
496:   return(0);
497: }

499: /*@
500:    ISGetSize - Returns the global length of an index set.

502:    Not Collective

504:    Input Parameter:
505: .  is - the index set

507:    Output Parameter:
508: .  size - the global size

510:    Level: beginner

512:    Concepts: size^of index set
513:    Concepts: index sets^size

515: @*/
516: PetscErrorCode  ISGetSize(IS is,PetscInt *size)
517: {

523:   (*is->ops->getsize)(is,size);
524:   return(0);
525: }

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

530:    Not Collective

532:    Input Parameter:
533: .  is - the index set

535:    Output Parameter:
536: .  size - the local size

538:    Level: beginner

540:    Concepts: size^of index set
541:    Concepts: local size^of index set
542:    Concepts: index sets^local size

544: @*/
545: PetscErrorCode  ISGetLocalSize(IS is,PetscInt *size)
546: {

552:   (*is->ops->getlocalsize)(is,size);
553:   return(0);
554: }

556: /*@C
557:    ISGetIndices - Returns a pointer to the indices.  The user should call
558:    ISRestoreIndices() after having looked at the indices.  The user should
559:    NOT change the indices.

561:    Not Collective

563:    Input Parameter:
564: .  is - the index set

566:    Output Parameter:
567: .  ptr - the location to put the pointer to the indices

569:    Fortran Note:
570:    This routine has two different interfaces from Fortran; the first is not recommend, it does not require Fortran 90
571: $    IS          is
572: $    integer     is_array(1)
573: $    PetscOffset i_is
574: $    int         ierr
575: $       call ISGetIndices(is,is_array,i_is,ierr)
576: $
577: $   Access first local entry in list
578: $      value = is_array(i_is + 1)
579: $
580: $      ...... other code
581: $       call ISRestoreIndices(is,is_array,i_is,ierr)
582:    The second Fortran interface is recommended.
583: $          use petscisdef
584: $          PetscInt, pointer :: array(:)
585: $          PetscErrorCode  ierr
586: $          IS       i
587: $          call ISGetIndicesF90(i,array,ierr)



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

594:    Level: intermediate

596:    Concepts: index sets^getting indices
597:    Concepts: indices of index set

599: .seealso: ISRestoreIndices(), ISGetIndicesF90()
600: @*/
601: PetscErrorCode  ISGetIndices(IS is,const PetscInt *ptr[])
602: {

608:   (*is->ops->getindices)(is,ptr);
609:   return(0);
610: }

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

615:    Not Collective

617:    Input Parameter:
618: .  is - the index set

620:    Output Parameter:
621: +   min - the minimum value
622: -   max - the maximum value

624:    Level: intermediate

626:    Notes:
627:     Empty index sets return min=PETSC_MAX_INT and max=PETSC_MIN_INT.
628:     In parallel, it returns the min and max of the local portion of the IS

630:    Concepts: index sets^getting indices
631:    Concepts: indices of index set

633: .seealso: ISGetIndices(), ISRestoreIndices(), ISGetIndicesF90()
634: @*/
635: PetscErrorCode  ISGetMinMax(IS is,PetscInt *min,PetscInt *max)
636: {
639:   if (min) *min = is->min;
640:   if (max) *max = is->max;
641:   return(0);
642: }

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

647:   Not Collective

649:   Input Parameter:
650: + is - the index set
651: - key - the search key

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

656:   Level: intermediate
657: @*/
658: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
659: {

663:   if (is->ops->locate) {
664:     (*is->ops->locate)(is,key,location);
665:   } else {
666:     PetscInt       numIdx;
667:     PetscBool      sorted;
668:     const PetscInt *idx;

670:     ISGetLocalSize(is,&numIdx);
671:     ISGetIndices(is,&idx);
672:     ISSorted(is,&sorted);
673:     if (sorted) {
674:       PetscFindInt(key,numIdx,idx,location);
675:     } else {
676:       PetscInt i;

678:       *location = -1;
679:       for (i = 0; i < numIdx; i++) {
680:         if (idx[i] == key) {
681:           *location = i;
682:           break;
683:         }
684:       }
685:     }
686:     ISRestoreIndices(is,&idx);
687:   }
688:   return(0);
689: }

691: /*@C
692:    ISRestoreIndices - Restores an index set to a usable state after a call
693:                       to ISGetIndices().

695:    Not Collective

697:    Input Parameters:
698: +  is - the index set
699: -  ptr - the pointer obtained by ISGetIndices()

701:    Fortran Note:
702:    This routine is used differently from Fortran
703: $    IS          is
704: $    integer     is_array(1)
705: $    PetscOffset i_is
706: $    int         ierr
707: $       call ISGetIndices(is,is_array,i_is,ierr)
708: $
709: $   Access first local entry in list
710: $      value = is_array(i_is + 1)
711: $
712: $      ...... other code
713: $       call ISRestoreIndices(is,is_array,i_is,ierr)

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

718:    Level: intermediate

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

723: .seealso: ISGetIndices(), ISRestoreIndicesF90()
724: @*/
725: PetscErrorCode  ISRestoreIndices(IS is,const PetscInt *ptr[])
726: {

732:   if (is->ops->restoreindices) {
733:     (*is->ops->restoreindices)(is,ptr);
734:   }
735:   return(0);
736: }

738: static PetscErrorCode ISGatherTotal_Private(IS is)
739: {
741:   PetscInt       i,n,N;
742:   const PetscInt *lindices;
743:   MPI_Comm       comm;
744:   PetscMPIInt    rank,size,*sizes = NULL,*offsets = NULL,nn;


749:   PetscObjectGetComm((PetscObject)is,&comm);
750:   MPI_Comm_size(comm,&size);
751:   MPI_Comm_rank(comm,&rank);
752:   ISGetLocalSize(is,&n);
753:   PetscMalloc2(size,&sizes,size,&offsets);

755:   PetscMPIIntCast(n,&nn);
756:   MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
757:   offsets[0] = 0;
758:   for (i=1; i<size; ++i) offsets[i] = offsets[i-1] + sizes[i-1];
759:   N = offsets[size-1] + sizes[size-1];

761:   PetscMalloc1(N,&(is->total));
762:   ISGetIndices(is,&lindices);
763:   MPI_Allgatherv((void*)lindices,nn,MPIU_INT,is->total,sizes,offsets,MPIU_INT,comm);
764:   ISRestoreIndices(is,&lindices);
765:   is->local_offset = offsets[rank];
766:   PetscFree2(sizes,offsets);
767:   return(0);
768: }

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

773:    Collective on IS

775:    Input Parameter:
776: .  is - the index set

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

782:    Level: intermediate

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

792:    Concepts: index sets^getting nonlocal indices
793: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices(), ISGetSize()
794: @*/
795: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
796: {
798:   PetscMPIInt    size;

803:   MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
804:   if (size == 1) {
805:     (*is->ops->getindices)(is,indices);
806:   } else {
807:     if (!is->total) {
808:       ISGatherTotal_Private(is);
809:     }
810:     *indices = is->total;
811:   }
812:   return(0);
813: }

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

818:    Not Collective.

820:    Input Parameter:
821: +  is - the index set
822: -  indices - index array; must be the array obtained with ISGetTotalIndices()

824:    Level: intermediate

826:    Concepts: index sets^getting nonlocal indices
827:    Concepts: index sets^restoring nonlocal indices
828: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices()
829: @*/
830: PetscErrorCode  ISRestoreTotalIndices(IS is, const PetscInt *indices[])
831: {
833:   PetscMPIInt    size;

838:   MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
839:   if (size == 1) {
840:     (*is->ops->restoreindices)(is,indices);
841:   } else {
842:     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.");
843:   }
844:   return(0);
845: }
846: /*@C
847:    ISGetNonlocalIndices - Retrieve an array of indices from remote processors
848:                        in this communicator.

850:    Collective on IS

852:    Input Parameter:
853: .  is - the index set

855:    Output Parameter:
856: .  indices - indices with rank 0 indices first, and so on,  omitting
857:              the current rank.  Total number of indices is the difference
858:              total and local, obtained with ISGetSize() and ISGetLocalSize(),
859:              respectively.

861:    Level: intermediate

863:    Notes:
864:     restore the indices using ISRestoreNonlocalIndices().
865:           The same scalability considerations as those for ISGetTotalIndices
866:           apply here.

868:    Concepts: index sets^getting nonlocal indices
869: .seealso: ISGetTotalIndices(), ISRestoreNonlocalIndices(), ISGetSize(), ISGetLocalSize().
870: @*/
871: PetscErrorCode  ISGetNonlocalIndices(IS is, const PetscInt *indices[])
872: {
874:   PetscMPIInt    size;
875:   PetscInt       n, N;

880:   MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
881:   if (size == 1) *indices = NULL;
882:   else {
883:     if (!is->total) {
884:       ISGatherTotal_Private(is);
885:     }
886:     ISGetLocalSize(is,&n);
887:     ISGetSize(is,&N);
888:     PetscMalloc1(N-n, &(is->nonlocal));
889:     PetscMemcpy(is->nonlocal, is->total, sizeof(PetscInt)*is->local_offset);
890:     PetscMemcpy(is->nonlocal+is->local_offset, is->total+is->local_offset+n, sizeof(PetscInt)*(N - is->local_offset - n));
891:     *indices = is->nonlocal;
892:   }
893:   return(0);
894: }

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

899:    Not Collective.

901:    Input Parameter:
902: +  is - the index set
903: -  indices - index array; must be the array obtained with ISGetNonlocalIndices()

905:    Level: intermediate

907:    Concepts: index sets^getting nonlocal indices
908:    Concepts: index sets^restoring nonlocal indices
909: .seealso: ISGetTotalIndices(), ISGetNonlocalIndices(), ISRestoreTotalIndices()
910: @*/
911: PetscErrorCode  ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
912: {
916:   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.");
917:   return(0);
918: }

920: /*@
921:    ISGetNonlocalIS - Gather all nonlocal indices for this IS and present
922:                      them as another sequential index set.


925:    Collective on IS

927:    Input Parameter:
928: .  is - the index set

930:    Output Parameter:
931: .  complement - sequential IS with indices identical to the result of
932:                 ISGetNonlocalIndices()

934:    Level: intermediate

936:    Notes:
937:     complement represents the result of ISGetNonlocalIndices as an IS.
938:           Therefore scalability issues similar to ISGetNonlocalIndices apply.
939:           The resulting IS must be restored using ISRestoreNonlocalIS().

941:    Concepts: index sets^getting nonlocal indices
942: .seealso: ISGetNonlocalIndices(), ISRestoreNonlocalIndices(),  ISAllGather(), ISGetSize()
943: @*/
944: PetscErrorCode  ISGetNonlocalIS(IS is, IS *complement)
945: {

951:   /* Check if the complement exists already. */
952:   if (is->complement) {
953:     *complement = is->complement;
954:     PetscObjectReference((PetscObject)(is->complement));
955:   } else {
956:     PetscInt       N, n;
957:     const PetscInt *idx;
958:     ISGetSize(is, &N);
959:     ISGetLocalSize(is,&n);
960:     ISGetNonlocalIndices(is, &idx);
961:     ISCreateGeneral(PETSC_COMM_SELF, N-n,idx, PETSC_USE_POINTER, &(is->complement));
962:     PetscObjectReference((PetscObject)is->complement);
963:     *complement = is->complement;
964:   }
965:   return(0);
966: }


969: /*@
970:    ISRestoreNonlocalIS - Restore the IS obtained with ISGetNonlocalIS().

972:    Not collective.

974:    Input Parameter:
975: +  is         - the index set
976: -  complement - index set of is's nonlocal indices

978:    Level: intermediate


981:    Concepts: index sets^getting nonlocal indices
982:    Concepts: index sets^restoring nonlocal indices
983: .seealso: ISGetNonlocalIS(), ISGetNonlocalIndices(), ISRestoreNonlocalIndices()
984: @*/
985: PetscErrorCode  ISRestoreNonlocalIS(IS is, IS *complement)
986: {
988:   PetscInt       refcnt;

993:   if (*complement != is->complement) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Complement IS being restored was not obtained with ISGetNonlocalIS()");
994:   PetscObjectGetReference((PetscObject)(is->complement), &refcnt);
995:   if (refcnt <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Duplicate call to ISRestoreNonlocalIS() detected");
996:   PetscObjectDereference((PetscObject)(is->complement));
997:   return(0);
998: }

1000: /*@C
1001:    ISView - Displays an index set.

1003:    Collective on IS

1005:    Input Parameters:
1006: +  is - the index set
1007: -  viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.

1009:    Level: intermediate

1011: .seealso: PetscViewerASCIIOpen()
1012: @*/
1013: PetscErrorCode  ISView(IS is,PetscViewer viewer)
1014: {

1019:   if (!viewer) {
1020:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is),&viewer);
1021:   }

1025:   PetscObjectPrintClassNamePrefixType((PetscObject)is,viewer);
1026:   (*is->ops->view)(is,viewer);
1027:   return(0);
1028: }

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

1033:   Collective on PetscViewer

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

1039:   Level: intermediate

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

1046:   Concepts: index set^loading from file

1048: .seealso: PetscViewerBinaryOpen(), ISView(), MatLoad(), VecLoad()
1049: @*/
1050: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
1051: {
1052:   PetscBool      isbinary, ishdf5;

1058:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
1059:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
1060:   if (!isbinary && !ishdf5) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1061:   if (!((PetscObject)is)->type_name) {ISSetType(is, ISGENERAL);}
1062:   (*is->ops->load)(is, viewer);
1063:   return(0);
1064: }

1066: /*@
1067:    ISSort - Sorts the indices of an index set.

1069:    Collective on IS

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

1074:    Level: intermediate

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

1079: .seealso: ISSortRemoveDups(), ISSorted()
1080: @*/
1081: PetscErrorCode  ISSort(IS is)
1082: {

1087:   (*is->ops->sort)(is);
1088:   return(0);
1089: }

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

1094:   Collective on IS

1096:   Input Parameters:
1097: . is - the index set

1099:   Level: intermediate

1101:   Concepts: index sets^sorting
1102:   Concepts: sorting^index set

1104: .seealso: ISSort(), ISSorted()
1105: @*/
1106: PetscErrorCode ISSortRemoveDups(IS is)
1107: {

1112:   (*is->ops->sortremovedups)(is);
1113:   return(0);
1114: }

1116: /*@
1117:    ISToGeneral - Converts an IS object of any type to ISGENERAL type

1119:    Collective on IS

1121:    Input Parameters:
1122: .  is - the index set

1124:    Level: intermediate

1126:    Concepts: index sets^sorting
1127:    Concepts: sorting^index set

1129: .seealso: ISSorted()
1130: @*/
1131: PetscErrorCode  ISToGeneral(IS is)
1132: {

1137:   if (is->ops->togeneral) {
1138:     (*is->ops->togeneral)(is);
1139:   } else SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Not written for this type %s",((PetscObject)is)->type_name);
1140:   return(0);
1141: }

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

1146:    Collective on IS

1148:    Input Parameter:
1149: .  is - the index set

1151:    Output Parameter:
1152: .  flg - output flag, either PETSC_TRUE if the index set is sorted,
1153:          or PETSC_FALSE otherwise.

1155:    Notes:
1156:     For parallel IS objects this only indicates if the local part of the IS
1157:           is sorted. So some processors may return PETSC_TRUE while others may
1158:           return PETSC_FALSE.

1160:    Level: intermediate

1162: .seealso: ISSort(), ISSortRemoveDups()
1163: @*/
1164: PetscErrorCode  ISSorted(IS is,PetscBool  *flg)
1165: {

1171:   (*is->ops->sorted)(is,flg);
1172:   return(0);
1173: }

1175: /*@
1176:    ISDuplicate - Creates a duplicate copy of an index set.

1178:    Collective on IS

1180:    Input Parmeters:
1181: .  is - the index set

1183:    Output Parameters:
1184: .  isnew - the copy of the index set

1186:    Level: beginner

1188:    Concepts: index sets^duplicating

1190: .seealso: ISCreateGeneral(), ISCopy()
1191: @*/
1192: PetscErrorCode  ISDuplicate(IS is,IS *newIS)
1193: {

1199:   (*is->ops->duplicate)(is,newIS);
1200:   (*newIS)->isidentity = is->isidentity;
1201:   (*newIS)->isperm     = is->isperm;
1202:   return(0);
1203: }

1205: /*@
1206:    ISCopy - Copies an index set.

1208:    Collective on IS

1210:    Input Parmeters:
1211: .  is - the index set

1213:    Output Parameters:
1214: .  isy - the copy of the index set

1216:    Level: beginner

1218:    Concepts: index sets^copying

1220: .seealso: ISDuplicate()
1221: @*/
1222: PetscErrorCode  ISCopy(IS is,IS isy)
1223: {

1230:   if (is == isy) return(0);
1231:   (*is->ops->copy)(is,isy);
1232:   isy->isperm     = is->isperm;
1233:   isy->max        = is->max;
1234:   isy->min        = is->min;
1235:   isy->isidentity = is->isidentity;
1236:   return(0);
1237: }

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

1242:    Collective on IS and comm

1244:    Input Arguments:
1245: + is - index set
1246: . comm - communicator for new index set
1247: - mode - copy semantics, PETSC_USE_POINTER for no-copy if possible, otherwise PETSC_COPY_VALUES

1249:    Output Arguments:
1250: . newis - new IS on comm

1252:    Level: advanced

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

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

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

1262: .seealso: ISSplit()
1263: @*/
1264: PetscErrorCode  ISOnComm(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
1265: {
1267:   PetscMPIInt    match;

1272:   MPI_Comm_compare(PetscObjectComm((PetscObject)is),comm,&match);
1273:   if (mode != PETSC_COPY_VALUES && (match == MPI_IDENT || match == MPI_CONGRUENT)) {
1274:     PetscObjectReference((PetscObject)is);
1275:     *newis = is;
1276:   } else {
1277:     (*is->ops->oncomm)(is,comm,mode,newis);
1278:   }
1279:   return(0);
1280: }

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

1285:    Logicall Collective on IS

1287:    Input Arguments:
1288: + is - index set
1289: - bs - block size

1291:    Level: intermediate

1293: .seealso: ISGetBlockSize(), ISCreateBlock()
1294: @*/
1295: PetscErrorCode  ISSetBlockSize(IS is,PetscInt bs)
1296: {

1302:   if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_OUTOFRANGE,"Block size %D, must be positive",bs);
1303:   (*is->ops->setblocksize)(is,bs);
1304:   return(0);
1305: }

1307: /*@
1308:    ISGetBlockSize - Returns the number of elements in a block.

1310:    Not Collective

1312:    Input Parameter:
1313: .  is - the index set

1315:    Output Parameter:
1316: .  size - the number of elements in a block

1318:    Level: intermediate

1320:    Concepts: IS^block size
1321:    Concepts: index sets^block size

1323: .seealso: ISBlockGetSize(), ISGetSize(), ISCreateBlock(), ISSetBlockSize()
1324: @*/
1325: PetscErrorCode  ISGetBlockSize(IS is,PetscInt *size)
1326: {

1330:   PetscLayoutGetBlockSize(is->map, size);
1331:   return(0);
1332: }

1334: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[])
1335: {
1337:   PetscInt       len,i;
1338:   const PetscInt *ptr;

1341:   ISGetSize(is,&len);
1342:   ISGetIndices(is,&ptr);
1343:   for (i=0; i<len; i++) idx[i] = ptr[i];
1344:   ISRestoreIndices(is,&ptr);
1345:   return(0);
1346: }

1348: /*MC
1349:     ISGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1350:     The users should call ISRestoreIndicesF90() after having looked at the
1351:     indices.  The user should NOT change the indices.

1353:     Synopsis:
1354:     ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1356:     Not collective

1358:     Input Parameter:
1359: .   x - index set

1361:     Output Parameters:
1362: +   xx_v - the Fortran90 pointer to the array
1363: -   ierr - error code

1365:     Example of Usage:
1366: .vb
1367:     PetscInt, 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:     Level: intermediate

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

1378:   Concepts: index sets^getting indices in f90
1379:   Concepts: indices of index set in f90

1381: M*/

1383: /*MC
1384:     ISRestoreIndicesF90 - Restores an index set to a usable state after
1385:     a call to ISGetIndicesF90().

1387:     Synopsis:
1388:     ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1390:     Not collective

1392:     Input Parameters:
1393: .   x - index set
1394: .   xx_v - the Fortran90 pointer to the array

1396:     Output Parameter:
1397: .   ierr - error code


1400:     Example of Usage:
1401: .vb
1402:     PetscInt, pointer xx_v(:)
1403:     ....
1404:     call ISGetIndicesF90(x,xx_v,ierr)
1405:     a = xx_v(3)
1406:     call ISRestoreIndicesF90(x,xx_v,ierr)
1407: .ve

1409:     Level: intermediate

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

1413: M*/

1415: /*MC
1416:     ISBlockGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1417:     The users should call ISBlockRestoreIndicesF90() after having looked at the
1418:     indices.  The user should NOT change the indices.

1420:     Synopsis:
1421:     ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1423:     Not collective

1425:     Input Parameter:
1426: .   x - index set

1428:     Output Parameters:
1429: +   xx_v - the Fortran90 pointer to the array
1430: -   ierr - error code
1431:     Example of Usage:
1432: .vb
1433:     PetscInt, pointer xx_v(:)
1434:     ....
1435:     call ISBlockGetIndicesF90(x,xx_v,ierr)
1436:     a = xx_v(3)
1437:     call ISBlockRestoreIndicesF90(x,xx_v,ierr)
1438: .ve

1440:     Level: intermediate

1442: .seealso:  ISBlockRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices(),
1443:            ISRestoreIndices()

1445:   Concepts: index sets^getting block indices in f90
1446:   Concepts: indices of index set in f90
1447:   Concepts: block^ indices of index set in f90

1449: M*/

1451: /*MC
1452:     ISBlockRestoreIndicesF90 - Restores an index set to a usable state after
1453:     a call to ISBlockGetIndicesF90().

1455:     Synopsis:
1456:     ISBlockRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)

1458:     Not Collective

1460:     Input Parameters:
1461: +   x - index set
1462: -   xx_v - the Fortran90 pointer to the array

1464:     Output Parameter:
1465: .   ierr - error code

1467:     Example of Usage:
1468: .vb
1469:     PetscInt, pointer xx_v(:)
1470:     ....
1471:     call ISBlockGetIndicesF90(x,xx_v,ierr)
1472:     a = xx_v(3)
1473:     call ISBlockRestoreIndicesF90(x,xx_v,ierr)
1474: .ve

1476:     Notes:
1477:     Not yet supported for all F90 compilers

1479:     Level: intermediate

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

1483: M*/