Actual source code: index.c
petsc-3.12.5 2020-03-29
2: /*
3: Defines the abstract operations on index sets, i.e. the public interface.
4: */
5: #include <petsc/private/isimpl.h>
6: #include <petscviewer.h>
7: #include <petscsf.h>
9: /* Logging support */
10: PetscClassId IS_CLASSID;
12: /*@
13: ISRenumber - Renumbers an index set (with multiplicities) in a contiguous way.
15: Collective on IS
17: Input Parmeters:
18: + subset - the index set
19: - subset_mult - the multiplcity of each entry in subset (optional, can be NULL)
21: Output Parameters:
22: + N - the maximum entry of the new IS
23: - subset_n - the new IS
25: Level: intermediate
27: .seealso:
28: @*/
29: PetscErrorCode ISRenumber(IS subset, IS subset_mult, PetscInt *N, IS *subset_n)
30: {
31: PetscSF sf;
32: PetscLayout map;
33: const PetscInt *idxs;
34: PetscInt *leaf_data,*root_data,*gidxs;
35: PetscInt N_n,n,i,lbounds[2],gbounds[2],Nl;
36: PetscInt n_n,nlocals,start,first_index;
37: PetscMPIInt commsize;
38: PetscBool first_found;
43: if (subset_mult) {
45: }
46: if (!N && !subset_n) return(0);
47: ISGetLocalSize(subset,&n);
48: if (subset_mult) {
49: ISGetLocalSize(subset_mult,&i);
50: if (i != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Local subset and multiplicity sizes don't match! %d != %d",n,i);
51: }
52: /* create workspace layout for computing global indices of subset */
53: ISGetIndices(subset,&idxs);
54: lbounds[0] = lbounds[1] = 0;
55: for (i=0;i<n;i++) {
56: if (idxs[i] < lbounds[0]) lbounds[0] = idxs[i];
57: else if (idxs[i] > lbounds[1]) lbounds[1] = idxs[i];
58: }
59: lbounds[0] = -lbounds[0];
60: MPIU_Allreduce(lbounds,gbounds,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)subset));
61: gbounds[0] = -gbounds[0];
62: N_n = gbounds[1] - gbounds[0] + 1;
64: PetscLayoutCreate(PetscObjectComm((PetscObject)subset),&map);
65: PetscLayoutSetBlockSize(map,1);
66: PetscLayoutSetSize(map,N_n);
67: PetscLayoutSetUp(map);
68: PetscLayoutGetLocalSize(map,&Nl);
70: /* create sf : leaf_data == multiplicity of indexes, root data == global index in layout */
71: PetscMalloc2(n,&leaf_data,Nl,&root_data);
72: if (subset_mult) {
73: const PetscInt* idxs_mult;
75: ISGetIndices(subset_mult,&idxs_mult);
76: PetscArraycpy(leaf_data,idxs_mult,n);
77: ISRestoreIndices(subset_mult,&idxs_mult);
78: } else {
79: for (i=0;i<n;i++) leaf_data[i] = 1;
80: }
81: /* local size of new subset */
82: n_n = 0;
83: for (i=0;i<n;i++) n_n += leaf_data[i];
85: /* global indexes in layout */
86: PetscMalloc1(n_n,&gidxs); /* allocating possibly extra space in gidxs which will be used later */
87: for (i=0;i<n;i++) gidxs[i] = idxs[i] - gbounds[0];
88: ISRestoreIndices(subset,&idxs);
89: PetscSFCreate(PetscObjectComm((PetscObject)subset),&sf);
90: PetscSFSetGraphLayout(sf,map,n,NULL,PETSC_COPY_VALUES,gidxs);
91: PetscLayoutDestroy(&map);
93: /* reduce from leaves to roots */
94: PetscArrayzero(root_data,Nl);
95: PetscSFReduceBegin(sf,MPIU_INT,leaf_data,root_data,MPI_MAX);
96: PetscSFReduceEnd(sf,MPIU_INT,leaf_data,root_data,MPI_MAX);
98: /* count indexes in local part of layout */
99: nlocals = 0;
100: first_index = -1;
101: first_found = PETSC_FALSE;
102: for (i=0;i<Nl;i++) {
103: if (!first_found && root_data[i]) {
104: first_found = PETSC_TRUE;
105: first_index = i;
106: }
107: nlocals += root_data[i];
108: }
110: /* cumulative of number of indexes and size of subset without holes */
111: #if defined(PETSC_HAVE_MPI_EXSCAN)
112: start = 0;
113: MPI_Exscan(&nlocals,&start,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)subset));
114: #else
115: MPI_Scan(&nlocals,&start,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)subset));
116: start = start-nlocals;
117: #endif
119: if (N) { /* compute total size of new subset if requested */
120: *N = start + nlocals;
121: MPI_Comm_size(PetscObjectComm((PetscObject)subset),&commsize);
122: MPI_Bcast(N,1,MPIU_INT,commsize-1,PetscObjectComm((PetscObject)subset));
123: }
125: if (!subset_n) {
126: PetscFree(gidxs);
127: PetscSFDestroy(&sf);
128: PetscFree2(leaf_data,root_data);
129: return(0);
130: }
132: /* adapt root data with cumulative */
133: if (first_found) {
134: PetscInt old_index;
136: root_data[first_index] += start;
137: old_index = first_index;
138: for (i=first_index+1;i<Nl;i++) {
139: if (root_data[i]) {
140: root_data[i] += root_data[old_index];
141: old_index = i;
142: }
143: }
144: }
146: /* from roots to leaves */
147: PetscSFBcastBegin(sf,MPIU_INT,root_data,leaf_data);
148: PetscSFBcastEnd(sf,MPIU_INT,root_data,leaf_data);
149: PetscSFDestroy(&sf);
151: /* create new IS with global indexes without holes */
152: if (subset_mult) {
153: const PetscInt* idxs_mult;
154: PetscInt cum;
156: cum = 0;
157: ISGetIndices(subset_mult,&idxs_mult);
158: for (i=0;i<n;i++) {
159: PetscInt j;
160: for (j=0;j<idxs_mult[i];j++) gidxs[cum++] = leaf_data[i] - idxs_mult[i] + j;
161: }
162: ISRestoreIndices(subset_mult,&idxs_mult);
163: } else {
164: for (i=0;i<n;i++) {
165: gidxs[i] = leaf_data[i]-1;
166: }
167: }
168: ISCreateGeneral(PetscObjectComm((PetscObject)subset),n_n,gidxs,PETSC_OWN_POINTER,subset_n);
169: PetscFree2(leaf_data,root_data);
170: return(0);
171: }
174: /*@
175: ISCreateSubIS - Create a sub index set from a global index set selecting some components.
177: Collective on IS
179: Input Parmeters:
180: + is - the index set
181: - comps - which components we will extract from is
183: Output Parameters:
184: . subis - the new sub index set
186: Level: intermediate
188: Example usage:
189: We have an index set (is) living on 3 processes with the following values:
190: | 4 9 0 | 2 6 7 | 10 11 1|
191: and another index set (comps) used to indicate which components of is we want to take,
192: | 7 5 | 1 2 | 0 4|
193: The output index set (subis) should look like:
194: | 11 7 | 9 0 | 4 6|
196: .seealso: VecGetSubVector(), MatCreateSubMatrix()
197: @*/
198: PetscErrorCode ISCreateSubIS(IS is,IS comps,IS *subis)
199: {
200: PetscSF sf;
201: const PetscInt *is_indices,*comps_indices;
202: PetscInt *subis_indices,nroots,nleaves,*mine,i,owner,lidx;
203: PetscSFNode *remote;
204: PetscErrorCode ierr;
205: MPI_Comm comm;
212: PetscObjectGetComm((PetscObject)is, &comm);
213: ISGetLocalSize(comps,&nleaves);
214: ISGetLocalSize(is,&nroots);
215: PetscMalloc1(nleaves,&remote);
216: PetscMalloc1(nleaves,&mine);
217: ISGetIndices(comps,&comps_indices);
218: /*
219: * Construct a PetscSF in which "is" data serves as roots and "subis" is leaves.
220: * Root data are sent to leaves using PetscSFBcast().
221: * */
222: for (i=0; i<nleaves; i++) {
223: mine[i] = i;
224: /* Connect a remote root with the current leaf. The value on the remote root
225: * will be received by the current local leaf.
226: * */
227: owner = -1;
228: lidx = -1;
229: PetscLayoutFindOwnerIndex(is->map,comps_indices[i],&owner, &lidx);
230: remote[i].rank = owner;
231: remote[i].index = lidx;
232: }
233: ISRestoreIndices(comps,&comps_indices);
234: PetscSFCreate(comm,&sf);
235: PetscSFSetFromOptions(sf);\
236: PetscSFSetGraph(sf,nroots,nleaves,mine,PETSC_OWN_POINTER,remote,PETSC_OWN_POINTER);
238: PetscMalloc1(nleaves,&subis_indices);
239: ISGetIndices(is, &is_indices);
240: PetscSFBcastBegin(sf,MPIU_INT,is_indices,subis_indices);
241: PetscSFBcastEnd(sf,MPIU_INT,is_indices,subis_indices);
242: ISRestoreIndices(is,&is_indices);
243: PetscSFDestroy(&sf);
244: ISCreateGeneral(comm,nleaves,subis_indices,PETSC_OWN_POINTER,subis);
245: return(0);
246: }
249: /*@
250: ISIdentity - Determines whether index set is the identity mapping.
252: Collective on IS
254: Input Parmeters:
255: . is - the index set
257: Output Parameters:
258: . ident - PETSC_TRUE if an identity, else PETSC_FALSE
260: Level: intermediate
263: .seealso: ISSetIdentity()
264: @*/
265: PetscErrorCode ISIdentity(IS is,PetscBool *ident)
266: {
272: *ident = is->isidentity;
273: if (*ident) return(0);
274: if (is->ops->identity) {
275: (*is->ops->identity)(is,ident);
276: }
277: return(0);
278: }
280: /*@
281: ISSetIdentity - Informs the index set that it is an identity.
283: Logically Collective on IS
285: Input Parmeters:
286: . is - the index set
288: Level: intermediate
291: .seealso: ISIdentity()
292: @*/
293: PetscErrorCode ISSetIdentity(IS is)
294: {
299: is->isidentity = PETSC_TRUE;
300: ISSetPermutation(is);
301: return(0);
302: }
304: /*@
305: ISContiguousLocal - Locates an index set with contiguous range within a global range, if possible
307: Not Collective
309: Input Parmeters:
310: + is - the index set
311: . gstart - global start
312: - gend - global end
314: Output Parameters:
315: + start - start of contiguous block, as an offset from gstart
316: - contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
318: Level: developer
320: .seealso: ISGetLocalSize(), VecGetOwnershipRange()
321: @*/
322: PetscErrorCode ISContiguousLocal(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
323: {
330: if (is->ops->contiguous) {
331: (*is->ops->contiguous)(is,gstart,gend,start,contig);
332: } else {
333: *start = -1;
334: *contig = PETSC_FALSE;
335: }
336: return(0);
337: }
339: /*@
340: ISPermutation - PETSC_TRUE or PETSC_FALSE depending on whether the
341: index set has been declared to be a permutation.
343: Logically Collective on IS
345: Input Parmeters:
346: . is - the index set
348: Output Parameters:
349: . perm - PETSC_TRUE if a permutation, else PETSC_FALSE
351: Level: intermediate
354: .seealso: ISSetPermutation()
355: @*/
356: PetscErrorCode ISPermutation(IS is,PetscBool *perm)
357: {
361: *perm = (PetscBool) is->isperm;
362: return(0);
363: }
365: /*@
366: ISSetPermutation - Informs the index set that it is a permutation.
368: Logically Collective on IS
370: Input Parmeters:
371: . is - the index set
373: Level: intermediate
376: The debug version of the libraries (./configure --with-debugging=1) checks if the
377: index set is actually a permutation. The optimized version just believes you.
379: .seealso: ISPermutation()
380: @*/
381: PetscErrorCode ISSetPermutation(IS is)
382: {
385: #if defined(PETSC_USE_DEBUG)
386: {
387: PetscMPIInt size;
390: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
391: if (size == 1) {
392: PetscInt i,n,*idx;
393: const PetscInt *iidx;
395: ISGetSize(is,&n);
396: PetscMalloc1(n,&idx);
397: ISGetIndices(is,&iidx);
398: PetscArraycpy(idx,iidx,n);
399: PetscSortInt(n,idx);
400: for (i=0; i<n; i++) {
401: if (idx[i] != i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Index set is not a permutation");
402: }
403: PetscFree(idx);
404: ISRestoreIndices(is,&iidx);
405: }
406: }
407: #endif
408: is->isperm = PETSC_TRUE;
409: return(0);
410: }
412: /*@
413: ISDestroy - Destroys an index set.
415: Collective on IS
417: Input Parameters:
418: . is - the index set
420: Level: beginner
422: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlocked()
423: @*/
424: PetscErrorCode ISDestroy(IS *is)
425: {
429: if (!*is) return(0);
431: if (--((PetscObject)(*is))->refct > 0) {*is = 0; return(0);}
432: if ((*is)->complement) {
433: PetscInt refcnt;
434: PetscObjectGetReference((PetscObject)((*is)->complement), &refcnt);
435: if (refcnt > 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Nonlocal IS has not been restored");
436: ISDestroy(&(*is)->complement);
437: }
438: if ((*is)->ops->destroy) {
439: (*(*is)->ops->destroy)(*is);
440: }
441: PetscLayoutDestroy(&(*is)->map);
442: /* Destroy local representations of offproc data. */
443: PetscFree((*is)->total);
444: PetscFree((*is)->nonlocal);
445: PetscHeaderDestroy(is);
446: return(0);
447: }
449: /*@
450: ISInvertPermutation - Creates a new permutation that is the inverse of
451: a given permutation.
453: Collective on IS
455: Input Parameter:
456: + is - the index set
457: - nlocal - number of indices on this processor in result (ignored for 1 proccessor) or
458: use PETSC_DECIDE
460: Output Parameter:
461: . isout - the inverse permutation
463: Level: intermediate
465: Notes:
466: For parallel index sets this does the complete parallel permutation, but the
467: code is not efficient for huge index sets (10,000,000 indices).
469: @*/
470: PetscErrorCode ISInvertPermutation(IS is,PetscInt nlocal,IS *isout)
471: {
477: if (!is->isperm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a permutation, must call ISSetPermutation() on the IS first");
478: if (is->isidentity) {
479: ISDuplicate(is,isout);
480: } else {
481: (*is->ops->invertpermutation)(is,nlocal,isout);
482: ISSetPermutation(*isout);
483: }
484: return(0);
485: }
487: /*@
488: ISGetSize - Returns the global length of an index set.
490: Not Collective
492: Input Parameter:
493: . is - the index set
495: Output Parameter:
496: . size - the global size
498: Level: beginner
501: @*/
502: PetscErrorCode ISGetSize(IS is,PetscInt *size)
503: {
507: *size = is->map->N;
508: return(0);
509: }
511: /*@
512: ISGetLocalSize - Returns the local (processor) length of an index set.
514: Not Collective
516: Input Parameter:
517: . is - the index set
519: Output Parameter:
520: . size - the local size
522: Level: beginner
524: @*/
525: PetscErrorCode ISGetLocalSize(IS is,PetscInt *size)
526: {
530: *size = is->map->n;
531: return(0);
532: }
534: /*@C
535: ISGetIndices - Returns a pointer to the indices. The user should call
536: ISRestoreIndices() after having looked at the indices. The user should
537: NOT change the indices.
539: Not Collective
541: Input Parameter:
542: . is - the index set
544: Output Parameter:
545: . ptr - the location to put the pointer to the indices
547: Fortran Note:
548: This routine has two different interfaces from Fortran; the first is not recommend, it does not require Fortran 90
549: $ IS is
550: $ integer is_array(1)
551: $ PetscOffset i_is
552: $ int ierr
553: $ call ISGetIndices(is,is_array,i_is,ierr)
554: $
555: $ Access first local entry in list
556: $ value = is_array(i_is + 1)
557: $
558: $ ...... other code
559: $ call ISRestoreIndices(is,is_array,i_is,ierr)
560: The second Fortran interface is recommended.
561: $ use petscisdef
562: $ PetscInt, pointer :: array(:)
563: $ PetscErrorCode ierr
564: $ IS i
565: $ call ISGetIndicesF90(i,array,ierr)
569: See the Fortran chapter of the users manual and
570: petsc/src/is/examples/[tutorials,tests] for details.
572: Level: intermediate
575: .seealso: ISRestoreIndices(), ISGetIndicesF90()
576: @*/
577: PetscErrorCode ISGetIndices(IS is,const PetscInt *ptr[])
578: {
584: (*is->ops->getindices)(is,ptr);
585: return(0);
586: }
588: /*@C
589: ISGetMinMax - Gets the minimum and maximum values in an IS
591: Not Collective
593: Input Parameter:
594: . is - the index set
596: Output Parameter:
597: + min - the minimum value
598: - max - the maximum value
600: Level: intermediate
602: Notes:
603: Empty index sets return min=PETSC_MAX_INT and max=PETSC_MIN_INT.
604: In parallel, it returns the min and max of the local portion of the IS
607: .seealso: ISGetIndices(), ISRestoreIndices(), ISGetIndicesF90()
608: @*/
609: PetscErrorCode ISGetMinMax(IS is,PetscInt *min,PetscInt *max)
610: {
613: if (min) *min = is->min;
614: if (max) *max = is->max;
615: return(0);
616: }
618: /*@
619: ISLocate - determine the location of an index within the local component of an index set
621: Not Collective
623: Input Parameter:
624: + is - the index set
625: - key - the search key
627: Output Parameter:
628: . location - if >= 0, a location within the index set that is equal to the key, otherwise the key is not in the index set
630: Level: intermediate
631: @*/
632: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
633: {
637: if (is->ops->locate) {
638: (*is->ops->locate)(is,key,location);
639: } else {
640: PetscInt numIdx;
641: PetscBool sorted;
642: const PetscInt *idx;
644: ISGetLocalSize(is,&numIdx);
645: ISGetIndices(is,&idx);
646: ISSorted(is,&sorted);
647: if (sorted) {
648: PetscFindInt(key,numIdx,idx,location);
649: } else {
650: PetscInt i;
652: *location = -1;
653: for (i = 0; i < numIdx; i++) {
654: if (idx[i] == key) {
655: *location = i;
656: break;
657: }
658: }
659: }
660: ISRestoreIndices(is,&idx);
661: }
662: return(0);
663: }
665: /*@C
666: ISRestoreIndices - Restores an index set to a usable state after a call
667: to ISGetIndices().
669: Not Collective
671: Input Parameters:
672: + is - the index set
673: - ptr - the pointer obtained by ISGetIndices()
675: Fortran Note:
676: This routine is used differently from Fortran
677: $ IS is
678: $ integer is_array(1)
679: $ PetscOffset i_is
680: $ int ierr
681: $ call ISGetIndices(is,is_array,i_is,ierr)
682: $
683: $ Access first local entry in list
684: $ value = is_array(i_is + 1)
685: $
686: $ ...... other code
687: $ call ISRestoreIndices(is,is_array,i_is,ierr)
689: See the Fortran chapter of the users manual and
690: petsc/src/is/examples/[tutorials,tests] for details.
692: Level: intermediate
694: Note:
695: This routine zeros out ptr. This is to prevent accidental us of the array after it has been restored.
697: .seealso: ISGetIndices(), ISRestoreIndicesF90()
698: @*/
699: PetscErrorCode ISRestoreIndices(IS is,const PetscInt *ptr[])
700: {
706: if (is->ops->restoreindices) {
707: (*is->ops->restoreindices)(is,ptr);
708: }
709: return(0);
710: }
712: static PetscErrorCode ISGatherTotal_Private(IS is)
713: {
715: PetscInt i,n,N;
716: const PetscInt *lindices;
717: MPI_Comm comm;
718: PetscMPIInt rank,size,*sizes = NULL,*offsets = NULL,nn;
723: PetscObjectGetComm((PetscObject)is,&comm);
724: MPI_Comm_size(comm,&size);
725: MPI_Comm_rank(comm,&rank);
726: ISGetLocalSize(is,&n);
727: PetscMalloc2(size,&sizes,size,&offsets);
729: PetscMPIIntCast(n,&nn);
730: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
731: offsets[0] = 0;
732: for (i=1; i<size; ++i) offsets[i] = offsets[i-1] + sizes[i-1];
733: N = offsets[size-1] + sizes[size-1];
735: PetscMalloc1(N,&(is->total));
736: ISGetIndices(is,&lindices);
737: MPI_Allgatherv((void*)lindices,nn,MPIU_INT,is->total,sizes,offsets,MPIU_INT,comm);
738: ISRestoreIndices(is,&lindices);
739: is->local_offset = offsets[rank];
740: PetscFree2(sizes,offsets);
741: return(0);
742: }
744: /*@C
745: ISGetTotalIndices - Retrieve an array containing all indices across the communicator.
747: Collective on IS
749: Input Parameter:
750: . is - the index set
752: Output Parameter:
753: . indices - total indices with rank 0 indices first, and so on; total array size is
754: the same as returned with ISGetSize().
756: Level: intermediate
758: Notes:
759: this is potentially nonscalable, but depends on the size of the total index set
760: and the size of the communicator. This may be feasible for index sets defined on
761: subcommunicators, such that the set size does not grow with PETSC_WORLD_COMM.
762: Note also that there is no way to tell where the local part of the indices starts
763: (use ISGetIndices() and ISGetNonlocalIndices() to retrieve just the local and just
764: the nonlocal part (complement), respectively).
766: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices(), ISGetSize()
767: @*/
768: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
769: {
771: PetscMPIInt size;
776: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
777: if (size == 1) {
778: (*is->ops->getindices)(is,indices);
779: } else {
780: if (!is->total) {
781: ISGatherTotal_Private(is);
782: }
783: *indices = is->total;
784: }
785: return(0);
786: }
788: /*@C
789: ISRestoreTotalIndices - Restore the index array obtained with ISGetTotalIndices().
791: Not Collective.
793: Input Parameter:
794: + is - the index set
795: - indices - index array; must be the array obtained with ISGetTotalIndices()
797: Level: intermediate
799: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices()
800: @*/
801: PetscErrorCode ISRestoreTotalIndices(IS is, const PetscInt *indices[])
802: {
804: PetscMPIInt size;
809: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
810: if (size == 1) {
811: (*is->ops->restoreindices)(is,indices);
812: } else {
813: if (is->total != *indices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index array pointer being restored does not point to the array obtained from the IS.");
814: }
815: return(0);
816: }
817: /*@C
818: ISGetNonlocalIndices - Retrieve an array of indices from remote processors
819: in this communicator.
821: Collective on IS
823: Input Parameter:
824: . is - the index set
826: Output Parameter:
827: . indices - indices with rank 0 indices first, and so on, omitting
828: the current rank. Total number of indices is the difference
829: total and local, obtained with ISGetSize() and ISGetLocalSize(),
830: respectively.
832: Level: intermediate
834: Notes:
835: restore the indices using ISRestoreNonlocalIndices().
836: The same scalability considerations as those for ISGetTotalIndices
837: apply here.
839: .seealso: ISGetTotalIndices(), ISRestoreNonlocalIndices(), ISGetSize(), ISGetLocalSize().
840: @*/
841: PetscErrorCode ISGetNonlocalIndices(IS is, const PetscInt *indices[])
842: {
844: PetscMPIInt size;
845: PetscInt n, N;
850: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
851: if (size == 1) *indices = NULL;
852: else {
853: if (!is->total) {
854: ISGatherTotal_Private(is);
855: }
856: ISGetLocalSize(is,&n);
857: ISGetSize(is,&N);
858: PetscMalloc1(N-n, &(is->nonlocal));
859: PetscArraycpy(is->nonlocal, is->total, is->local_offset);
860: PetscArraycpy(is->nonlocal+is->local_offset, is->total+is->local_offset+n,N - is->local_offset - n);
861: *indices = is->nonlocal;
862: }
863: return(0);
864: }
866: /*@C
867: ISRestoreTotalIndices - Restore the index array obtained with ISGetNonlocalIndices().
869: Not Collective.
871: Input Parameter:
872: + is - the index set
873: - indices - index array; must be the array obtained with ISGetNonlocalIndices()
875: Level: intermediate
877: .seealso: ISGetTotalIndices(), ISGetNonlocalIndices(), ISRestoreTotalIndices()
878: @*/
879: PetscErrorCode ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
880: {
884: if (is->nonlocal != *indices) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index array pointer being restored does not point to the array obtained from the IS.");
885: return(0);
886: }
888: /*@
889: ISGetNonlocalIS - Gather all nonlocal indices for this IS and present
890: them as another sequential index set.
893: Collective on IS
895: Input Parameter:
896: . is - the index set
898: Output Parameter:
899: . complement - sequential IS with indices identical to the result of
900: ISGetNonlocalIndices()
902: Level: intermediate
904: Notes:
905: complement represents the result of ISGetNonlocalIndices as an IS.
906: Therefore scalability issues similar to ISGetNonlocalIndices apply.
907: The resulting IS must be restored using ISRestoreNonlocalIS().
909: .seealso: ISGetNonlocalIndices(), ISRestoreNonlocalIndices(), ISAllGather(), ISGetSize()
910: @*/
911: PetscErrorCode ISGetNonlocalIS(IS is, IS *complement)
912: {
918: /* Check if the complement exists already. */
919: if (is->complement) {
920: *complement = is->complement;
921: PetscObjectReference((PetscObject)(is->complement));
922: } else {
923: PetscInt N, n;
924: const PetscInt *idx;
925: ISGetSize(is, &N);
926: ISGetLocalSize(is,&n);
927: ISGetNonlocalIndices(is, &idx);
928: ISCreateGeneral(PETSC_COMM_SELF, N-n,idx, PETSC_USE_POINTER, &(is->complement));
929: PetscObjectReference((PetscObject)is->complement);
930: *complement = is->complement;
931: }
932: return(0);
933: }
936: /*@
937: ISRestoreNonlocalIS - Restore the IS obtained with ISGetNonlocalIS().
939: Not collective.
941: Input Parameter:
942: + is - the index set
943: - complement - index set of is's nonlocal indices
945: Level: intermediate
948: .seealso: ISGetNonlocalIS(), ISGetNonlocalIndices(), ISRestoreNonlocalIndices()
949: @*/
950: PetscErrorCode ISRestoreNonlocalIS(IS is, IS *complement)
951: {
953: PetscInt refcnt;
958: if (*complement != is->complement) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Complement IS being restored was not obtained with ISGetNonlocalIS()");
959: PetscObjectGetReference((PetscObject)(is->complement), &refcnt);
960: if (refcnt <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Duplicate call to ISRestoreNonlocalIS() detected");
961: PetscObjectDereference((PetscObject)(is->complement));
962: return(0);
963: }
965: /*@C
966: ISView - Displays an index set.
968: Collective on IS
970: Input Parameters:
971: + is - the index set
972: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
974: Level: intermediate
976: .seealso: PetscViewerASCIIOpen()
977: @*/
978: PetscErrorCode ISView(IS is,PetscViewer viewer)
979: {
984: if (!viewer) {
985: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is),&viewer);
986: }
990: PetscObjectPrintClassNamePrefixType((PetscObject)is,viewer);
991: (*is->ops->view)(is,viewer);
992: return(0);
993: }
995: /*@
996: ISLoad - Loads a vector that has been stored in binary or HDF5 format with ISView().
998: Collective on PetscViewer
1000: Input Parameters:
1001: + is - the newly loaded vector, this needs to have been created with ISCreate() or some related function before a call to ISLoad().
1002: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or HDF5 file viewer, obtained from PetscViewerHDF5Open()
1004: Level: intermediate
1006: Notes:
1007: IF using HDF5, you must assign the IS the same name as was used in the IS
1008: that was stored in the file using PetscObjectSetName(). Otherwise you will
1009: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"
1011: .seealso: PetscViewerBinaryOpen(), ISView(), MatLoad(), VecLoad()
1012: @*/
1013: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
1014: {
1015: PetscBool isbinary, ishdf5;
1021: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
1022: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
1023: if (!isbinary && !ishdf5) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1024: if (!((PetscObject)is)->type_name) {ISSetType(is, ISGENERAL);}
1025: (*is->ops->load)(is, viewer);
1026: return(0);
1027: }
1029: /*@
1030: ISSort - Sorts the indices of an index set.
1032: Collective on IS
1034: Input Parameters:
1035: . is - the index set
1037: Level: intermediate
1040: .seealso: ISSortRemoveDups(), ISSorted()
1041: @*/
1042: PetscErrorCode ISSort(IS is)
1043: {
1048: (*is->ops->sort)(is);
1049: return(0);
1050: }
1052: /*@
1053: ISSortRemoveDups - Sorts the indices of an index set, removing duplicates.
1055: Collective on IS
1057: Input Parameters:
1058: . is - the index set
1060: Level: intermediate
1063: .seealso: ISSort(), ISSorted()
1064: @*/
1065: PetscErrorCode ISSortRemoveDups(IS is)
1066: {
1071: (*is->ops->sortremovedups)(is);
1072: return(0);
1073: }
1075: /*@
1076: ISToGeneral - Converts an IS object of any type to ISGENERAL type
1078: Collective on IS
1080: Input Parameters:
1081: . is - the index set
1083: Level: intermediate
1086: .seealso: ISSorted()
1087: @*/
1088: PetscErrorCode ISToGeneral(IS is)
1089: {
1094: if (is->ops->togeneral) {
1095: (*is->ops->togeneral)(is);
1096: } else SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Not written for this type %s",((PetscObject)is)->type_name);
1097: return(0);
1098: }
1100: /*@
1101: ISSorted - Checks the indices to determine whether they have been sorted.
1103: Collective on IS
1105: Input Parameter:
1106: . is - the index set
1108: Output Parameter:
1109: . flg - output flag, either PETSC_TRUE if the index set is sorted,
1110: or PETSC_FALSE otherwise.
1112: Notes:
1113: For parallel IS objects this only indicates if the local part of the IS
1114: is sorted. So some processors may return PETSC_TRUE while others may
1115: return PETSC_FALSE.
1117: Level: intermediate
1119: .seealso: ISSort(), ISSortRemoveDups()
1120: @*/
1121: PetscErrorCode ISSorted(IS is,PetscBool *flg)
1122: {
1128: (*is->ops->sorted)(is,flg);
1129: return(0);
1130: }
1132: /*@
1133: ISDuplicate - Creates a duplicate copy of an index set.
1135: Collective on IS
1137: Input Parmeters:
1138: . is - the index set
1140: Output Parameters:
1141: . isnew - the copy of the index set
1143: Level: beginner
1145: .seealso: ISCreateGeneral(), ISCopy()
1146: @*/
1147: PetscErrorCode ISDuplicate(IS is,IS *newIS)
1148: {
1154: (*is->ops->duplicate)(is,newIS);
1155: (*newIS)->isidentity = is->isidentity;
1156: (*newIS)->isperm = is->isperm;
1157: return(0);
1158: }
1160: /*@
1161: ISCopy - Copies an index set.
1163: Collective on IS
1165: Input Parmeters:
1166: . is - the index set
1168: Output Parameters:
1169: . isy - the copy of the index set
1171: Level: beginner
1173: .seealso: ISDuplicate()
1174: @*/
1175: PetscErrorCode ISCopy(IS is,IS isy)
1176: {
1183: if (is == isy) return(0);
1184: (*is->ops->copy)(is,isy);
1185: isy->isperm = is->isperm;
1186: isy->max = is->max;
1187: isy->min = is->min;
1188: isy->isidentity = is->isidentity;
1189: return(0);
1190: }
1192: /*@
1193: ISOnComm - Split a parallel IS on subcomms (usually self) or concatenate index sets on subcomms into a parallel index set
1195: Collective on IS
1197: Input Arguments:
1198: + is - index set
1199: . comm - communicator for new index set
1200: - mode - copy semantics, PETSC_USE_POINTER for no-copy if possible, otherwise PETSC_COPY_VALUES
1202: Output Arguments:
1203: . newis - new IS on comm
1205: Level: advanced
1207: Notes:
1208: It is usually desirable to create a parallel IS and look at the local part when necessary.
1210: This function is useful if serial ISs must be created independently, or to view many
1211: logically independent serial ISs.
1213: The input IS must have the same type on every process.
1215: .seealso: ISSplit()
1216: @*/
1217: PetscErrorCode ISOnComm(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
1218: {
1220: PetscMPIInt match;
1225: MPI_Comm_compare(PetscObjectComm((PetscObject)is),comm,&match);
1226: if (mode != PETSC_COPY_VALUES && (match == MPI_IDENT || match == MPI_CONGRUENT)) {
1227: PetscObjectReference((PetscObject)is);
1228: *newis = is;
1229: } else {
1230: (*is->ops->oncomm)(is,comm,mode,newis);
1231: }
1232: return(0);
1233: }
1235: /*@
1236: ISSetBlockSize - informs an index set that it has a given block size
1238: Logicall Collective on IS
1240: Input Arguments:
1241: + is - index set
1242: - bs - block size
1244: Level: intermediate
1246: .seealso: ISGetBlockSize(), ISCreateBlock()
1247: @*/
1248: PetscErrorCode ISSetBlockSize(IS is,PetscInt bs)
1249: {
1255: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_OUTOFRANGE,"Block size %D, must be positive",bs);
1256: (*is->ops->setblocksize)(is,bs);
1257: return(0);
1258: }
1260: /*@
1261: ISGetBlockSize - Returns the number of elements in a block.
1263: Not Collective
1265: Input Parameter:
1266: . is - the index set
1268: Output Parameter:
1269: . size - the number of elements in a block
1271: Level: intermediate
1274: .seealso: ISBlockGetSize(), ISGetSize(), ISCreateBlock(), ISSetBlockSize()
1275: @*/
1276: PetscErrorCode ISGetBlockSize(IS is,PetscInt *size)
1277: {
1281: PetscLayoutGetBlockSize(is->map, size);
1282: return(0);
1283: }
1285: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[])
1286: {
1288: PetscInt len,i;
1289: const PetscInt *ptr;
1292: ISGetSize(is,&len);
1293: ISGetIndices(is,&ptr);
1294: for (i=0; i<len; i++) idx[i] = ptr[i];
1295: ISRestoreIndices(is,&ptr);
1296: return(0);
1297: }
1299: /*MC
1300: ISGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1301: The users should call ISRestoreIndicesF90() after having looked at the
1302: indices. The user should NOT change the indices.
1304: Synopsis:
1305: ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1307: Not collective
1309: Input Parameter:
1310: . x - index set
1312: Output Parameters:
1313: + xx_v - the Fortran90 pointer to the array
1314: - ierr - error code
1316: Example of Usage:
1317: .vb
1318: PetscInt, pointer xx_v(:)
1319: ....
1320: call ISGetIndicesF90(x,xx_v,ierr)
1321: a = xx_v(3)
1322: call ISRestoreIndicesF90(x,xx_v,ierr)
1323: .ve
1325: Level: intermediate
1327: .seealso: ISRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices()
1330: M*/
1332: /*MC
1333: ISRestoreIndicesF90 - Restores an index set to a usable state after
1334: a call to ISGetIndicesF90().
1336: Synopsis:
1337: ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1339: Not collective
1341: Input Parameters:
1342: + x - index set
1343: - xx_v - the Fortran90 pointer to the array
1345: Output Parameter:
1346: . ierr - error code
1349: Example of Usage:
1350: .vb
1351: PetscInt, pointer xx_v(:)
1352: ....
1353: call ISGetIndicesF90(x,xx_v,ierr)
1354: a = xx_v(3)
1355: call ISRestoreIndicesF90(x,xx_v,ierr)
1356: .ve
1358: Level: intermediate
1360: .seealso: ISGetIndicesF90(), ISGetIndices(), ISRestoreIndices()
1362: M*/
1364: /*MC
1365: ISBlockGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1366: The users should call ISBlockRestoreIndicesF90() after having looked at the
1367: indices. The user should NOT change the indices.
1369: Synopsis:
1370: ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1372: Not collective
1374: Input Parameter:
1375: . x - index set
1377: Output Parameters:
1378: + xx_v - the Fortran90 pointer to the array
1379: - ierr - error code
1380: Example of Usage:
1381: .vb
1382: PetscInt, pointer xx_v(:)
1383: ....
1384: call ISBlockGetIndicesF90(x,xx_v,ierr)
1385: a = xx_v(3)
1386: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
1387: .ve
1389: Level: intermediate
1391: .seealso: ISBlockRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices(),
1392: ISRestoreIndices()
1394: M*/
1396: /*MC
1397: ISBlockRestoreIndicesF90 - Restores an index set to a usable state after
1398: a call to ISBlockGetIndicesF90().
1400: Synopsis:
1401: ISBlockRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1403: Not Collective
1405: Input Parameters:
1406: + x - index set
1407: - xx_v - the Fortran90 pointer to the array
1409: Output Parameter:
1410: . ierr - error code
1412: Example of Usage:
1413: .vb
1414: PetscInt, pointer xx_v(:)
1415: ....
1416: call ISBlockGetIndicesF90(x,xx_v,ierr)
1417: a = xx_v(3)
1418: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
1419: .ve
1421: Notes:
1422: Not yet supported for all F90 compilers
1424: Level: intermediate
1426: .seealso: ISBlockGetIndicesF90(), ISGetIndices(), ISRestoreIndices(), ISRestoreIndicesF90()
1428: M*/