Actual source code: index.c
petsc-3.9.4 2018-09-11
2: /*
3: Defines the abstract operations on index sets, i.e. the public interface.
4: */
5: #include <petsc/private/isimpl.h>
6: #include <petscviewer.h>
7: #include <petscsf.h>
9: /* Logging support */
10: PetscClassId IS_CLASSID;
12: /*@
13: ISRenumber - Renumbers an index set (with multiplicities) in a contiguous way.
15: Collective on IS
17: Input Parmeters:
18: . subset - the index set
19: . subset_mult - the multiplcity of each entry in subset (optional, can be NULL)
21: Output Parameters:
22: . N - the maximum entry of the new IS
23: . subset_n - the new IS
25: Level: intermediate
27: .seealso:
28: @*/
29: PetscErrorCode ISRenumber(IS subset, IS subset_mult, PetscInt *N, IS *subset_n)
30: {
31: PetscSF sf;
32: PetscLayout map;
33: const PetscInt *idxs;
34: PetscInt *leaf_data,*root_data,*gidxs;
35: PetscInt N_n,n,i,lbounds[2],gbounds[2],Nl;
36: PetscInt n_n,nlocals,start,first_index;
37: PetscMPIInt commsize;
38: PetscBool first_found;
42: ISGetLocalSize(subset,&n);
43: if (subset_mult) {
45: ISGetLocalSize(subset,&i);
46: if (i != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Local subset and multiplicity sizes don't match! %d != %d",n,i);
47: }
48: /* create workspace layout for computing global indices of subset */
49: ISGetIndices(subset,&idxs);
50: lbounds[0] = lbounds[1] = 0;
51: for (i=0;i<n;i++) {
52: if (idxs[i] < lbounds[0]) lbounds[0] = idxs[i];
53: else if (idxs[i] > lbounds[1]) lbounds[1] = idxs[i];
54: }
55: lbounds[0] = -lbounds[0];
56: MPIU_Allreduce(lbounds,gbounds,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)subset));
57: gbounds[0] = -gbounds[0];
58: N_n= gbounds[1] - gbounds[0] + 1;
59: PetscLayoutCreate(PetscObjectComm((PetscObject)subset),&map);
60: PetscLayoutSetBlockSize(map,1);
61: PetscLayoutSetSize(map,N_n);
62: PetscLayoutSetUp(map);
63: PetscLayoutGetLocalSize(map,&Nl);
65: /* create sf : leaf_data == multiplicity of indexes, root data == global index in layout */
66: PetscMalloc2(n,&leaf_data,Nl,&root_data);
67: if (subset_mult) {
68: const PetscInt* idxs_mult;
70: ISGetIndices(subset_mult,&idxs_mult);
71: PetscMemcpy(leaf_data,idxs_mult,n*sizeof(PetscInt));
72: ISRestoreIndices(subset_mult,&idxs_mult);
73: } else {
74: for (i=0;i<n;i++) leaf_data[i] = 1;
75: }
76: /* local size of new subset */
77: n_n = 0;
78: for (i=0;i<n;i++) n_n += leaf_data[i];
80: /* global indexes in layout */
81: PetscMalloc1(n_n,&gidxs); /* allocating possibly extra space in gidxs which will be used later */
82: for (i=0;i<n;i++) gidxs[i] = idxs[i] - gbounds[0];
83: ISRestoreIndices(subset,&idxs);
84: PetscSFCreate(PetscObjectComm((PetscObject)subset),&sf);
85: PetscSFSetGraphLayout(sf,map,n,NULL,PETSC_COPY_VALUES,gidxs);
86: PetscLayoutDestroy(&map);
88: /* reduce from leaves to roots */
89: PetscMemzero(root_data,Nl*sizeof(PetscInt));
90: PetscSFReduceBegin(sf,MPIU_INT,leaf_data,root_data,MPI_MAX);
91: PetscSFReduceEnd(sf,MPIU_INT,leaf_data,root_data,MPI_MAX);
93: /* count indexes in local part of layout */
94: nlocals = 0;
95: first_index = -1;
96: first_found = PETSC_FALSE;
97: for (i=0;i<Nl;i++) {
98: if (!first_found && root_data[i]) {
99: first_found = PETSC_TRUE;
100: first_index = i;
101: }
102: nlocals += root_data[i];
103: }
105: /* cumulative of number of indexes and size of subset without holes */
106: #if defined(PETSC_HAVE_MPI_EXSCAN)
107: start = 0;
108: MPI_Exscan(&nlocals,&start,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)subset));
109: #else
110: MPI_Scan(&nlocals,&start,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)subset));
111: start = start-nlocals;
112: #endif
114: if (N) { /* compute total size of new subset if requested */
115: *N = start + nlocals;
116: MPI_Comm_size(PetscObjectComm((PetscObject)subset),&commsize);
117: MPI_Bcast(N,1,MPIU_INT,commsize-1,PetscObjectComm((PetscObject)subset));
118: }
120: /* adapt root data with cumulative */
121: if (first_found) {
122: PetscInt old_index;
124: root_data[first_index] += start;
125: old_index = first_index;
126: for (i=first_index+1;i<Nl;i++) {
127: if (root_data[i]) {
128: root_data[i] += root_data[old_index];
129: old_index = i;
130: }
131: }
132: }
134: /* from roots to leaves */
135: PetscSFBcastBegin(sf,MPIU_INT,root_data,leaf_data);
136: PetscSFBcastEnd(sf,MPIU_INT,root_data,leaf_data);
137: PetscSFDestroy(&sf);
139: /* create new IS with global indexes without holes */
140: if (subset_mult) {
141: const PetscInt* idxs_mult;
142: PetscInt cum;
144: cum = 0;
145: ISGetIndices(subset_mult,&idxs_mult);
146: for (i=0;i<n;i++) {
147: PetscInt j;
148: for (j=0;j<idxs_mult[i];j++) gidxs[cum++] = leaf_data[i] - idxs_mult[i] + j;
149: }
150: ISRestoreIndices(subset_mult,&idxs_mult);
151: } else {
152: for (i=0;i<n;i++) {
153: gidxs[i] = leaf_data[i]-1;
154: }
155: }
156: ISCreateGeneral(PetscObjectComm((PetscObject)subset),n_n,gidxs,PETSC_OWN_POINTER,subset_n);
157: PetscFree2(leaf_data,root_data);
158: return(0);
159: }
161: /*@
162: ISIdentity - Determines whether index set is the identity mapping.
164: Collective on IS
166: Input Parmeters:
167: . is - the index set
169: Output Parameters:
170: . ident - PETSC_TRUE if an identity, else PETSC_FALSE
172: Level: intermediate
174: Concepts: identity mapping
175: Concepts: index sets^is identity
177: .seealso: ISSetIdentity()
178: @*/
179: PetscErrorCode ISIdentity(IS is,PetscBool *ident)
180: {
186: *ident = is->isidentity;
187: if (*ident) return(0);
188: if (is->ops->identity) {
189: (*is->ops->identity)(is,ident);
190: }
191: return(0);
192: }
194: /*@
195: ISSetIdentity - Informs the index set that it is an identity.
197: Logically Collective on IS
199: Input Parmeters:
200: . is - the index set
202: Level: intermediate
204: Concepts: identity mapping
205: Concepts: index sets^is identity
207: .seealso: ISIdentity()
208: @*/
209: PetscErrorCode ISSetIdentity(IS is)
210: {
215: is->isidentity = PETSC_TRUE;
216: ISSetPermutation(is);
217: return(0);
218: }
220: /*@
221: ISContiguousLocal - Locates an index set with contiguous range within a global range, if possible
223: Not Collective
225: Input Parmeters:
226: + is - the index set
227: . gstart - global start
228: . gend - global end
230: Output Parameters:
231: + start - start of contiguous block, as an offset from gstart
232: - contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
234: Level: developer
236: Concepts: index sets^is contiguous
238: .seealso: ISGetLocalSize(), VecGetOwnershipRange()
239: @*/
240: PetscErrorCode ISContiguousLocal(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
241: {
248: if (is->ops->contiguous) {
249: (*is->ops->contiguous)(is,gstart,gend,start,contig);
250: } else {
251: *start = -1;
252: *contig = PETSC_FALSE;
253: }
254: return(0);
255: }
257: /*@
258: ISPermutation - PETSC_TRUE or PETSC_FALSE depending on whether the
259: index set has been declared to be a permutation.
261: Logically Collective on IS
263: Input Parmeters:
264: . is - the index set
266: Output Parameters:
267: . perm - PETSC_TRUE if a permutation, else PETSC_FALSE
269: Level: intermediate
271: Concepts: permutation
272: Concepts: index sets^is permutation
274: .seealso: ISSetPermutation()
275: @*/
276: PetscErrorCode ISPermutation(IS is,PetscBool *perm)
277: {
281: *perm = (PetscBool) is->isperm;
282: return(0);
283: }
285: /*@
286: ISSetPermutation - Informs the index set that it is a permutation.
288: Logically Collective on IS
290: Input Parmeters:
291: . is - the index set
293: Level: intermediate
295: Concepts: permutation
296: Concepts: index sets^permutation
298: The debug version of the libraries (./configure --with-debugging=1) checks if the
299: index set is actually a permutation. The optimized version just believes you.
301: .seealso: ISPermutation()
302: @*/
303: PetscErrorCode ISSetPermutation(IS is)
304: {
307: #if defined(PETSC_USE_DEBUG)
308: {
309: PetscMPIInt size;
312: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
313: if (size == 1) {
314: PetscInt i,n,*idx;
315: const PetscInt *iidx;
317: ISGetSize(is,&n);
318: PetscMalloc1(n,&idx);
319: ISGetIndices(is,&iidx);
320: PetscMemcpy(idx,iidx,n*sizeof(PetscInt));
321: PetscSortInt(n,idx);
322: for (i=0; i<n; i++) {
323: if (idx[i] != i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Index set is not a permutation");
324: }
325: PetscFree(idx);
326: ISRestoreIndices(is,&iidx);
327: }
328: }
329: #endif
330: is->isperm = PETSC_TRUE;
331: return(0);
332: }
334: /*@
335: ISDestroy - Destroys an index set.
337: Collective on IS
339: Input Parameters:
340: . is - the index set
342: Level: beginner
344: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlocked()
345: @*/
346: PetscErrorCode ISDestroy(IS *is)
347: {
351: if (!*is) return(0);
353: if (--((PetscObject)(*is))->refct > 0) {*is = 0; return(0);}
354: if ((*is)->complement) {
355: PetscInt refcnt;
356: PetscObjectGetReference((PetscObject)((*is)->complement), &refcnt);
357: if (refcnt > 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Nonlocal IS has not been restored");
358: ISDestroy(&(*is)->complement);
359: }
360: if ((*is)->ops->destroy) {
361: (*(*is)->ops->destroy)(*is);
362: }
363: PetscLayoutDestroy(&(*is)->map);
364: /* Destroy local representations of offproc data. */
365: PetscFree((*is)->total);
366: PetscFree((*is)->nonlocal);
367: PetscHeaderDestroy(is);
368: return(0);
369: }
371: /*@
372: ISInvertPermutation - Creates a new permutation that is the inverse of
373: a given permutation.
375: Collective on IS
377: Input Parameter:
378: + is - the index set
379: - nlocal - number of indices on this processor in result (ignored for 1 proccessor) or
380: use PETSC_DECIDE
382: Output Parameter:
383: . isout - the inverse permutation
385: Level: intermediate
387: Notes: For parallel index sets this does the complete parallel permutation, but the
388: code is not efficient for huge index sets (10,000,000 indices).
390: Concepts: inverse permutation
391: Concepts: permutation^inverse
392: Concepts: index sets^inverting
393: @*/
394: PetscErrorCode ISInvertPermutation(IS is,PetscInt nlocal,IS *isout)
395: {
401: if (!is->isperm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a permutation, must call ISSetPermutation() on the IS first");
402: if (is->isidentity) {
403: ISDuplicate(is,isout);
404: } else {
405: (*is->ops->invertpermutation)(is,nlocal,isout);
406: ISSetPermutation(*isout);
407: }
408: return(0);
409: }
411: /*@
412: ISGetSize - Returns the global length of an index set.
414: Not Collective
416: Input Parameter:
417: . is - the index set
419: Output Parameter:
420: . size - the global size
422: Level: beginner
424: Concepts: size^of index set
425: Concepts: index sets^size
427: @*/
428: PetscErrorCode ISGetSize(IS is,PetscInt *size)
429: {
435: (*is->ops->getsize)(is,size);
436: return(0);
437: }
439: /*@
440: ISGetLocalSize - Returns the local (processor) length of an index set.
442: Not Collective
444: Input Parameter:
445: . is - the index set
447: Output Parameter:
448: . size - the local size
450: Level: beginner
452: Concepts: size^of index set
453: Concepts: local size^of index set
454: Concepts: index sets^local size
456: @*/
457: PetscErrorCode ISGetLocalSize(IS is,PetscInt *size)
458: {
464: (*is->ops->getlocalsize)(is,size);
465: return(0);
466: }
468: /*@C
469: ISGetIndices - Returns a pointer to the indices. The user should call
470: ISRestoreIndices() after having looked at the indices. The user should
471: NOT change the indices.
473: Not Collective
475: Input Parameter:
476: . is - the index set
478: Output Parameter:
479: . ptr - the location to put the pointer to the indices
481: Fortran Note:
482: This routine has two different interfaces from Fortran; the first is not recommend, it does not require Fortran 90
483: $ IS is
484: $ integer is_array(1)
485: $ PetscOffset i_is
486: $ int ierr
487: $ call ISGetIndices(is,is_array,i_is,ierr)
488: $
489: $ Access first local entry in list
490: $ value = is_array(i_is + 1)
491: $
492: $ ...... other code
493: $ call ISRestoreIndices(is,is_array,i_is,ierr)
494: The second Fortran interface is recommended.
495: $ use petscisdef
496: $ PetscInt, pointer :: array(:)
497: $ PetscErrorCode ierr
498: $ IS i
499: $ call ISGetIndicesF90(i,array,ierr)
503: See the Fortran chapter of the users manual and
504: petsc/src/is/examples/[tutorials,tests] for details.
506: Level: intermediate
508: Concepts: index sets^getting indices
509: Concepts: indices of index set
511: .seealso: ISRestoreIndices(), ISGetIndicesF90()
512: @*/
513: PetscErrorCode ISGetIndices(IS is,const PetscInt *ptr[])
514: {
520: (*is->ops->getindices)(is,ptr);
521: return(0);
522: }
524: /*@C
525: ISGetMinMax - Gets the minimum and maximum values in an IS
527: Not Collective
529: Input Parameter:
530: . is - the index set
532: Output Parameter:
533: + min - the minimum value
534: - max - the maximum value
536: Level: intermediate
538: Notes:
539: Empty index sets return min=PETSC_MAX_INT and max=PETSC_MIN_INT.
540: In parallel, it returns the min and max of the local portion of the IS
542: Concepts: index sets^getting indices
543: Concepts: indices of index set
545: .seealso: ISGetIndices(), ISRestoreIndices(), ISGetIndicesF90()
546: @*/
547: PetscErrorCode ISGetMinMax(IS is,PetscInt *min,PetscInt *max)
548: {
551: if (min) *min = is->min;
552: if (max) *max = is->max;
553: return(0);
554: }
556: /*@
557: ISLocate - determine the location of an index within the local component of an index set
559: Not Collective
561: Input Parameter:
562: + is - the index set
563: - key - the search key
565: Output Parameter:
566: . location - if >= 0, a location within the index set that is equal to the key, otherwise the key is not in the index set
568: Level: intermediate
569: @*/
570: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
571: {
575: if (is->ops->locate) {
576: (*is->ops->locate)(is,key,location);
577: } else {
578: PetscInt numIdx;
579: PetscBool sorted;
580: const PetscInt *idx;
582: ISGetLocalSize(is,&numIdx);
583: ISGetIndices(is,&idx);
584: ISSorted(is,&sorted);
585: if (sorted) {
586: PetscFindInt(key,numIdx,idx,location);
587: } else {
588: PetscInt i;
590: *location = -1;
591: for (i = 0; i < numIdx; i++) {
592: if (idx[i] == key) {
593: *location = i;
594: break;
595: }
596: }
597: }
598: ISRestoreIndices(is,&idx);
599: }
600: return(0);
601: }
603: /*@C
604: ISRestoreIndices - Restores an index set to a usable state after a call
605: to ISGetIndices().
607: Not Collective
609: Input Parameters:
610: + is - the index set
611: - ptr - the pointer obtained by ISGetIndices()
613: Fortran Note:
614: This routine is used differently from Fortran
615: $ IS is
616: $ integer is_array(1)
617: $ PetscOffset i_is
618: $ int ierr
619: $ call ISGetIndices(is,is_array,i_is,ierr)
620: $
621: $ Access first local entry in list
622: $ value = is_array(i_is + 1)
623: $
624: $ ...... other code
625: $ call ISRestoreIndices(is,is_array,i_is,ierr)
627: See the Fortran chapter of the users manual and
628: petsc/src/is/examples/[tutorials,tests] for details.
630: Level: intermediate
632: Note:
633: This routine zeros out ptr. This is to prevent accidental us of the array after it has been restored.
635: .seealso: ISGetIndices(), ISRestoreIndicesF90()
636: @*/
637: PetscErrorCode ISRestoreIndices(IS is,const PetscInt *ptr[])
638: {
644: if (is->ops->restoreindices) {
645: (*is->ops->restoreindices)(is,ptr);
646: }
647: return(0);
648: }
650: static PetscErrorCode ISGatherTotal_Private(IS is)
651: {
653: PetscInt i,n,N;
654: const PetscInt *lindices;
655: MPI_Comm comm;
656: PetscMPIInt rank,size,*sizes = NULL,*offsets = NULL,nn;
661: PetscObjectGetComm((PetscObject)is,&comm);
662: MPI_Comm_size(comm,&size);
663: MPI_Comm_rank(comm,&rank);
664: ISGetLocalSize(is,&n);
665: PetscMalloc2(size,&sizes,size,&offsets);
667: PetscMPIIntCast(n,&nn);
668: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
669: offsets[0] = 0;
670: for (i=1; i<size; ++i) offsets[i] = offsets[i-1] + sizes[i-1];
671: N = offsets[size-1] + sizes[size-1];
673: PetscMalloc1(N,&(is->total));
674: ISGetIndices(is,&lindices);
675: MPI_Allgatherv((void*)lindices,nn,MPIU_INT,is->total,sizes,offsets,MPIU_INT,comm);
676: ISRestoreIndices(is,&lindices);
677: is->local_offset = offsets[rank];
678: PetscFree2(sizes,offsets);
679: return(0);
680: }
682: /*@C
683: ISGetTotalIndices - Retrieve an array containing all indices across the communicator.
685: Collective on IS
687: Input Parameter:
688: . is - the index set
690: Output Parameter:
691: . indices - total indices with rank 0 indices first, and so on; total array size is
692: the same as returned with ISGetSize().
694: Level: intermediate
696: Notes: this is potentially nonscalable, but depends on the size of the total index set
697: and the size of the communicator. This may be feasible for index sets defined on
698: subcommunicators, such that the set size does not grow with PETSC_WORLD_COMM.
699: Note also that there is no way to tell where the local part of the indices starts
700: (use ISGetIndices() and ISGetNonlocalIndices() to retrieve just the local and just
701: the nonlocal part (complement), respectively).
703: Concepts: index sets^getting nonlocal indices
704: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices(), ISGetSize()
705: @*/
706: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
707: {
709: PetscMPIInt size;
714: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
715: if (size == 1) {
716: (*is->ops->getindices)(is,indices);
717: } else {
718: if (!is->total) {
719: ISGatherTotal_Private(is);
720: }
721: *indices = is->total;
722: }
723: return(0);
724: }
726: /*@C
727: ISRestoreTotalIndices - Restore the index array obtained with ISGetTotalIndices().
729: Not Collective.
731: Input Parameter:
732: + is - the index set
733: - indices - index array; must be the array obtained with ISGetTotalIndices()
735: Level: intermediate
737: Concepts: index sets^getting nonlocal indices
738: Concepts: index sets^restoring nonlocal indices
739: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices()
740: @*/
741: PetscErrorCode ISRestoreTotalIndices(IS is, const PetscInt *indices[])
742: {
744: PetscMPIInt size;
749: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
750: if (size == 1) {
751: (*is->ops->restoreindices)(is,indices);
752: } else {
753: if (is->total != *indices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index array pointer being restored does not point to the array obtained from the IS.");
754: }
755: return(0);
756: }
757: /*@C
758: ISGetNonlocalIndices - Retrieve an array of indices from remote processors
759: in this communicator.
761: Collective on IS
763: Input Parameter:
764: . is - the index set
766: Output Parameter:
767: . indices - indices with rank 0 indices first, and so on, omitting
768: the current rank. Total number of indices is the difference
769: total and local, obtained with ISGetSize() and ISGetLocalSize(),
770: respectively.
772: Level: intermediate
774: Notes: restore the indices using ISRestoreNonlocalIndices().
775: The same scalability considerations as those for ISGetTotalIndices
776: apply here.
778: Concepts: index sets^getting nonlocal indices
779: .seealso: ISGetTotalIndices(), ISRestoreNonlocalIndices(), ISGetSize(), ISGetLocalSize().
780: @*/
781: PetscErrorCode ISGetNonlocalIndices(IS is, const PetscInt *indices[])
782: {
784: PetscMPIInt size;
785: PetscInt n, N;
790: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
791: if (size == 1) *indices = NULL;
792: else {
793: if (!is->total) {
794: ISGatherTotal_Private(is);
795: }
796: ISGetLocalSize(is,&n);
797: ISGetSize(is,&N);
798: PetscMalloc1(N-n, &(is->nonlocal));
799: PetscMemcpy(is->nonlocal, is->total, sizeof(PetscInt)*is->local_offset);
800: PetscMemcpy(is->nonlocal+is->local_offset, is->total+is->local_offset+n, sizeof(PetscInt)*(N - is->local_offset - n));
801: *indices = is->nonlocal;
802: }
803: return(0);
804: }
806: /*@C
807: ISRestoreTotalIndices - Restore the index array obtained with ISGetNonlocalIndices().
809: Not Collective.
811: Input Parameter:
812: + is - the index set
813: - indices - index array; must be the array obtained with ISGetNonlocalIndices()
815: Level: intermediate
817: Concepts: index sets^getting nonlocal indices
818: Concepts: index sets^restoring nonlocal indices
819: .seealso: ISGetTotalIndices(), ISGetNonlocalIndices(), ISRestoreTotalIndices()
820: @*/
821: PetscErrorCode ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
822: {
826: if (is->nonlocal != *indices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index array pointer being restored does not point to the array obtained from the IS.");
827: return(0);
828: }
830: /*@
831: ISGetNonlocalIS - Gather all nonlocal indices for this IS and present
832: them as another sequential index set.
835: Collective on IS
837: Input Parameter:
838: . is - the index set
840: Output Parameter:
841: . complement - sequential IS with indices identical to the result of
842: ISGetNonlocalIndices()
844: Level: intermediate
846: Notes: complement represents the result of ISGetNonlocalIndices as an IS.
847: Therefore scalability issues similar to ISGetNonlocalIndices apply.
848: The resulting IS must be restored using ISRestoreNonlocalIS().
850: Concepts: index sets^getting nonlocal indices
851: .seealso: ISGetNonlocalIndices(), ISRestoreNonlocalIndices(), ISAllGather(), ISGetSize()
852: @*/
853: PetscErrorCode ISGetNonlocalIS(IS is, IS *complement)
854: {
860: /* Check if the complement exists already. */
861: if (is->complement) {
862: *complement = is->complement;
863: PetscObjectReference((PetscObject)(is->complement));
864: } else {
865: PetscInt N, n;
866: const PetscInt *idx;
867: ISGetSize(is, &N);
868: ISGetLocalSize(is,&n);
869: ISGetNonlocalIndices(is, &idx);
870: ISCreateGeneral(PETSC_COMM_SELF, N-n,idx, PETSC_USE_POINTER, &(is->complement));
871: PetscObjectReference((PetscObject)is->complement);
872: *complement = is->complement;
873: }
874: return(0);
875: }
878: /*@
879: ISRestoreNonlocalIS - Restore the IS obtained with ISGetNonlocalIS().
881: Not collective.
883: Input Parameter:
884: + is - the index set
885: - complement - index set of is's nonlocal indices
887: Level: intermediate
890: Concepts: index sets^getting nonlocal indices
891: Concepts: index sets^restoring nonlocal indices
892: .seealso: ISGetNonlocalIS(), ISGetNonlocalIndices(), ISRestoreNonlocalIndices()
893: @*/
894: PetscErrorCode ISRestoreNonlocalIS(IS is, IS *complement)
895: {
897: PetscInt refcnt;
902: if (*complement != is->complement) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Complement IS being restored was not obtained with ISGetNonlocalIS()");
903: PetscObjectGetReference((PetscObject)(is->complement), &refcnt);
904: if (refcnt <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Duplicate call to ISRestoreNonlocalIS() detected");
905: PetscObjectDereference((PetscObject)(is->complement));
906: return(0);
907: }
909: /*@C
910: ISView - Displays an index set.
912: Collective on IS
914: Input Parameters:
915: + is - the index set
916: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
918: Level: intermediate
920: .seealso: PetscViewerASCIIOpen()
921: @*/
922: PetscErrorCode ISView(IS is,PetscViewer viewer)
923: {
928: if (!viewer) {
929: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is),&viewer);
930: }
934: PetscObjectPrintClassNamePrefixType((PetscObject)is,viewer);
935: (*is->ops->view)(is,viewer);
936: return(0);
937: }
939: /*@
940: ISLoad - Loads a vector that has been stored in binary or HDF5 format with ISView().
942: Collective on PetscViewer
944: Input Parameters:
945: + is - the newly loaded vector, this needs to have been created with ISCreate() or some related function before a call to ISLoad().
946: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or HDF5 file viewer, obtained from PetscViewerHDF5Open()
948: Level: intermediate
950: Notes:
951: IF using HDF5, you must assign the IS the same name as was used in the IS
952: that was stored in the file using PetscObjectSetName(). Otherwise you will
953: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"
955: Concepts: index set^loading from file
957: .seealso: PetscViewerBinaryOpen(), ISView(), MatLoad(), VecLoad()
958: @*/
959: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
960: {
961: PetscBool isbinary, ishdf5;
967: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
968: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
969: if (!isbinary && !ishdf5) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
970: if (!((PetscObject)is)->type_name) {ISSetType(is, ISGENERAL);}
971: (*is->ops->load)(is, viewer);
972: return(0);
973: }
975: /*@
976: ISSort - Sorts the indices of an index set.
978: Collective on IS
980: Input Parameters:
981: . is - the index set
983: Level: intermediate
985: Concepts: index sets^sorting
986: Concepts: sorting^index set
988: .seealso: ISSortRemoveDups(), ISSorted()
989: @*/
990: PetscErrorCode ISSort(IS is)
991: {
996: (*is->ops->sort)(is);
997: return(0);
998: }
1000: /*@
1001: ISSortRemoveDups - Sorts the indices of an index set, removing duplicates.
1003: Collective on IS
1005: Input Parameters:
1006: . is - the index set
1008: Level: intermediate
1010: Concepts: index sets^sorting
1011: Concepts: sorting^index set
1013: .seealso: ISSort(), ISSorted()
1014: @*/
1015: PetscErrorCode ISSortRemoveDups(IS is)
1016: {
1021: (*is->ops->sortremovedups)(is);
1022: return(0);
1023: }
1025: /*@
1026: ISToGeneral - Converts an IS object of any type to ISGENERAL type
1028: Collective on IS
1030: Input Parameters:
1031: . is - the index set
1033: Level: intermediate
1035: Concepts: index sets^sorting
1036: Concepts: sorting^index set
1038: .seealso: ISSorted()
1039: @*/
1040: PetscErrorCode ISToGeneral(IS is)
1041: {
1046: if (is->ops->togeneral) {
1047: (*is->ops->togeneral)(is);
1048: } else SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Not written for this type %s",((PetscObject)is)->type_name);
1049: return(0);
1050: }
1052: /*@
1053: ISSorted - Checks the indices to determine whether they have been sorted.
1055: Collective on IS
1057: Input Parameter:
1058: . is - the index set
1060: Output Parameter:
1061: . flg - output flag, either PETSC_TRUE if the index set is sorted,
1062: or PETSC_FALSE otherwise.
1064: Notes: For parallel IS objects this only indicates if the local part of the IS
1065: is sorted. So some processors may return PETSC_TRUE while others may
1066: return PETSC_FALSE.
1068: Level: intermediate
1070: .seealso: ISSort(), ISSortRemoveDups()
1071: @*/
1072: PetscErrorCode ISSorted(IS is,PetscBool *flg)
1073: {
1079: (*is->ops->sorted)(is,flg);
1080: return(0);
1081: }
1083: /*@
1084: ISDuplicate - Creates a duplicate copy of an index set.
1086: Collective on IS
1088: Input Parmeters:
1089: . is - the index set
1091: Output Parameters:
1092: . isnew - the copy of the index set
1094: Level: beginner
1096: Concepts: index sets^duplicating
1098: .seealso: ISCreateGeneral(), ISCopy()
1099: @*/
1100: PetscErrorCode ISDuplicate(IS is,IS *newIS)
1101: {
1107: (*is->ops->duplicate)(is,newIS);
1108: (*newIS)->isidentity = is->isidentity;
1109: (*newIS)->isperm = is->isperm;
1110: return(0);
1111: }
1113: /*@
1114: ISCopy - Copies an index set.
1116: Collective on IS
1118: Input Parmeters:
1119: . is - the index set
1121: Output Parameters:
1122: . isy - the copy of the index set
1124: Level: beginner
1126: Concepts: index sets^copying
1128: .seealso: ISDuplicate()
1129: @*/
1130: PetscErrorCode ISCopy(IS is,IS isy)
1131: {
1138: if (is == isy) return(0);
1139: (*is->ops->copy)(is,isy);
1140: isy->isperm = is->isperm;
1141: isy->max = is->max;
1142: isy->min = is->min;
1143: isy->isidentity = is->isidentity;
1144: return(0);
1145: }
1147: /*@
1148: ISOnComm - Split a parallel IS on subcomms (usually self) or concatenate index sets on subcomms into a parallel index set
1150: Collective on IS and comm
1152: Input Arguments:
1153: + is - index set
1154: . comm - communicator for new index set
1155: - mode - copy semantics, PETSC_USE_POINTER for no-copy if possible, otherwise PETSC_COPY_VALUES
1157: Output Arguments:
1158: . newis - new IS on comm
1160: Level: advanced
1162: Notes:
1163: It is usually desirable to create a parallel IS and look at the local part when necessary.
1165: This function is useful if serial ISs must be created independently, or to view many
1166: logically independent serial ISs.
1168: The input IS must have the same type on every process.
1170: .seealso: ISSplit()
1171: @*/
1172: PetscErrorCode ISOnComm(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
1173: {
1175: PetscMPIInt match;
1180: MPI_Comm_compare(PetscObjectComm((PetscObject)is),comm,&match);
1181: if (mode != PETSC_COPY_VALUES && (match == MPI_IDENT || match == MPI_CONGRUENT)) {
1182: PetscObjectReference((PetscObject)is);
1183: *newis = is;
1184: } else {
1185: (*is->ops->oncomm)(is,comm,mode,newis);
1186: }
1187: return(0);
1188: }
1190: /*@
1191: ISSetBlockSize - informs an index set that it has a given block size
1193: Logicall Collective on IS
1195: Input Arguments:
1196: + is - index set
1197: - bs - block size
1199: Level: intermediate
1201: .seealso: ISGetBlockSize(), ISCreateBlock()
1202: @*/
1203: PetscErrorCode ISSetBlockSize(IS is,PetscInt bs)
1204: {
1210: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_OUTOFRANGE,"Block size %D, must be positive",bs);
1211: (*is->ops->setblocksize)(is,bs);
1212: return(0);
1213: }
1215: /*@
1216: ISGetBlockSize - Returns the number of elements in a block.
1218: Not Collective
1220: Input Parameter:
1221: . is - the index set
1223: Output Parameter:
1224: . size - the number of elements in a block
1226: Level: intermediate
1228: Concepts: IS^block size
1229: Concepts: index sets^block size
1231: .seealso: ISBlockGetSize(), ISGetSize(), ISCreateBlock(), ISSetBlockSize()
1232: @*/
1233: PetscErrorCode ISGetBlockSize(IS is,PetscInt *size)
1234: {
1238: PetscLayoutGetBlockSize(is->map, size);
1239: return(0);
1240: }
1242: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[])
1243: {
1245: PetscInt len,i;
1246: const PetscInt *ptr;
1249: ISGetSize(is,&len);
1250: ISGetIndices(is,&ptr);
1251: for (i=0; i<len; i++) idx[i] = ptr[i];
1252: ISRestoreIndices(is,&ptr);
1253: return(0);
1254: }
1256: /*MC
1257: ISGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1258: The users should call ISRestoreIndicesF90() after having looked at the
1259: indices. The user should NOT change the indices.
1261: Synopsis:
1262: ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1264: Not collective
1266: Input Parameter:
1267: . x - index set
1269: Output Parameters:
1270: + xx_v - the Fortran90 pointer to the array
1271: - ierr - error code
1273: Example of Usage:
1274: .vb
1275: PetscScalar, pointer xx_v(:)
1276: ....
1277: call ISGetIndicesF90(x,xx_v,ierr)
1278: a = xx_v(3)
1279: call ISRestoreIndicesF90(x,xx_v,ierr)
1280: .ve
1282: Notes:
1283: Not yet supported for all F90 compilers.
1285: Level: intermediate
1287: .seealso: ISRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices()
1289: Concepts: index sets^getting indices in f90
1290: Concepts: indices of index set in f90
1292: M*/
1294: /*MC
1295: ISRestoreIndicesF90 - Restores an index set to a usable state after
1296: a call to ISGetIndicesF90().
1298: Synopsis:
1299: ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1301: Not collective
1303: Input Parameters:
1304: . x - index set
1305: . xx_v - the Fortran90 pointer to the array
1307: Output Parameter:
1308: . ierr - error code
1311: Example of Usage:
1312: .vb
1313: PetscScalar, pointer xx_v(:)
1314: ....
1315: call ISGetIndicesF90(x,xx_v,ierr)
1316: a = xx_v(3)
1317: call ISRestoreIndicesF90(x,xx_v,ierr)
1318: .ve
1320: Notes:
1321: Not yet supported for all F90 compilers.
1323: Level: intermediate
1325: .seealso: ISGetIndicesF90(), ISGetIndices(), ISRestoreIndices()
1327: M*/
1329: /*MC
1330: ISBlockGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1331: The users should call ISBlockRestoreIndicesF90() after having looked at the
1332: indices. The user should NOT change the indices.
1334: Synopsis:
1335: ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1337: Not collective
1339: Input Parameter:
1340: . x - index set
1342: Output Parameters:
1343: + xx_v - the Fortran90 pointer to the array
1344: - ierr - error code
1345: Example of Usage:
1346: .vb
1347: PetscScalar, pointer xx_v(:)
1348: ....
1349: call ISBlockGetIndicesF90(x,xx_v,ierr)
1350: a = xx_v(3)
1351: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
1352: .ve
1354: Notes:
1355: Not yet supported for all F90 compilers
1357: Level: intermediate
1359: .seealso: ISBlockRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices(),
1360: ISRestoreIndices()
1362: Concepts: index sets^getting block indices in f90
1363: Concepts: indices of index set in f90
1364: Concepts: block^ indices of index set in f90
1366: M*/
1368: /*MC
1369: ISBlockRestoreIndicesF90 - Restores an index set to a usable state after
1370: a call to ISBlockGetIndicesF90().
1372: Synopsis:
1373: ISBlockRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1375: Not Collective
1377: Input Parameters:
1378: + x - index set
1379: - xx_v - the Fortran90 pointer to the array
1381: Output Parameter:
1382: . ierr - error code
1384: Example of Usage:
1385: .vb
1386: PetscScalar, pointer xx_v(:)
1387: ....
1388: call ISBlockGetIndicesF90(x,xx_v,ierr)
1389: a = xx_v(3)
1390: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
1391: .ve
1393: Notes:
1394: Not yet supported for all F90 compilers
1396: Level: intermediate
1398: .seealso: ISBlockGetIndicesF90(), ISGetIndices(), ISRestoreIndices(), ISRestoreIndicesF90()
1400: M*/