Actual source code: index.c
1: /*
2: Defines the abstract operations on index sets, i.e. the public interface.
3: */
4: #include <petsc/private/isimpl.h>
5: #include <petscviewer.h>
6: #include <petscsf.h>
8: /* Logging support */
9: PetscClassId IS_CLASSID;
10: /* TODO: Much more events are missing! */
11: PetscLogEvent IS_View;
12: PetscLogEvent IS_Load;
14: /*@
15: ISRenumber - Renumbers an index set (with multiplicities) in a contiguous way.
17: Collective on IS
19: Input Parmeters:
20: + subset - the index set
21: - subset_mult - the multiplcity of each entry in subset (optional, can be NULL)
23: Output Parameters:
24: + N - the maximum entry of the new IS
25: - subset_n - the new IS
27: Level: intermediate
29: .seealso:
30: @*/
31: PetscErrorCode ISRenumber(IS subset, IS subset_mult, PetscInt *N, IS *subset_n)
32: {
33: PetscSF sf;
34: PetscLayout map;
35: const PetscInt *idxs;
36: PetscInt *leaf_data,*root_data,*gidxs;
37: PetscInt N_n,n,i,lbounds[2],gbounds[2],Nl;
38: PetscInt n_n,nlocals,start,first_index;
39: PetscMPIInt commsize;
40: PetscBool first_found;
45: if (subset_mult) {
47: }
48: if (!N && !subset_n) return(0);
49: ISGetLocalSize(subset,&n);
50: if (subset_mult) {
51: ISGetLocalSize(subset_mult,&i);
52: if (i != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Local subset and multiplicity sizes don't match! %d != %d",n,i);
53: }
54: /* create workspace layout for computing global indices of subset */
55: ISGetIndices(subset,&idxs);
56: lbounds[0] = lbounds[1] = 0;
57: for (i=0;i<n;i++) {
58: if (idxs[i] < lbounds[0]) lbounds[0] = idxs[i];
59: else if (idxs[i] > lbounds[1]) lbounds[1] = idxs[i];
60: }
61: lbounds[0] = -lbounds[0];
62: MPIU_Allreduce(lbounds,gbounds,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)subset));
63: gbounds[0] = -gbounds[0];
64: N_n = gbounds[1] - gbounds[0] + 1;
66: PetscLayoutCreate(PetscObjectComm((PetscObject)subset),&map);
67: PetscLayoutSetBlockSize(map,1);
68: PetscLayoutSetSize(map,N_n);
69: PetscLayoutSetUp(map);
70: PetscLayoutGetLocalSize(map,&Nl);
72: /* create sf : leaf_data == multiplicity of indexes, root data == global index in layout */
73: PetscMalloc2(n,&leaf_data,Nl,&root_data);
74: if (subset_mult) {
75: const PetscInt* idxs_mult;
77: ISGetIndices(subset_mult,&idxs_mult);
78: PetscArraycpy(leaf_data,idxs_mult,n);
79: ISRestoreIndices(subset_mult,&idxs_mult);
80: } else {
81: for (i=0;i<n;i++) leaf_data[i] = 1;
82: }
83: /* local size of new subset */
84: n_n = 0;
85: for (i=0;i<n;i++) n_n += leaf_data[i];
87: /* global indexes in layout */
88: PetscMalloc1(n_n,&gidxs); /* allocating possibly extra space in gidxs which will be used later */
89: for (i=0;i<n;i++) gidxs[i] = idxs[i] - gbounds[0];
90: ISRestoreIndices(subset,&idxs);
91: PetscSFCreate(PetscObjectComm((PetscObject)subset),&sf);
92: PetscSFSetGraphLayout(sf,map,n,NULL,PETSC_COPY_VALUES,gidxs);
93: PetscLayoutDestroy(&map);
95: /* reduce from leaves to roots */
96: PetscArrayzero(root_data,Nl);
97: PetscSFReduceBegin(sf,MPIU_INT,leaf_data,root_data,MPI_MAX);
98: PetscSFReduceEnd(sf,MPIU_INT,leaf_data,root_data,MPI_MAX);
100: /* count indexes in local part of layout */
101: nlocals = 0;
102: first_index = -1;
103: first_found = PETSC_FALSE;
104: for (i=0;i<Nl;i++) {
105: if (!first_found && root_data[i]) {
106: first_found = PETSC_TRUE;
107: first_index = i;
108: }
109: nlocals += root_data[i];
110: }
112: /* cumulative of number of indexes and size of subset without holes */
113: #if defined(PETSC_HAVE_MPI_EXSCAN)
114: start = 0;
115: MPI_Exscan(&nlocals,&start,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)subset));
116: #else
117: MPI_Scan(&nlocals,&start,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)subset));
118: start = start-nlocals;
119: #endif
121: if (N) { /* compute total size of new subset if requested */
122: *N = start + nlocals;
123: MPI_Comm_size(PetscObjectComm((PetscObject)subset),&commsize);
124: MPI_Bcast(N,1,MPIU_INT,commsize-1,PetscObjectComm((PetscObject)subset));
125: }
127: if (!subset_n) {
128: PetscFree(gidxs);
129: PetscSFDestroy(&sf);
130: PetscFree2(leaf_data,root_data);
131: return(0);
132: }
134: /* adapt root data with cumulative */
135: if (first_found) {
136: PetscInt old_index;
138: root_data[first_index] += start;
139: old_index = first_index;
140: for (i=first_index+1;i<Nl;i++) {
141: if (root_data[i]) {
142: root_data[i] += root_data[old_index];
143: old_index = i;
144: }
145: }
146: }
148: /* from roots to leaves */
149: PetscSFBcastBegin(sf,MPIU_INT,root_data,leaf_data,MPI_REPLACE);
150: PetscSFBcastEnd(sf,MPIU_INT,root_data,leaf_data,MPI_REPLACE);
151: PetscSFDestroy(&sf);
153: /* create new IS with global indexes without holes */
154: if (subset_mult) {
155: const PetscInt* idxs_mult;
156: PetscInt cum;
158: cum = 0;
159: ISGetIndices(subset_mult,&idxs_mult);
160: for (i=0;i<n;i++) {
161: PetscInt j;
162: for (j=0;j<idxs_mult[i];j++) gidxs[cum++] = leaf_data[i] - idxs_mult[i] + j;
163: }
164: ISRestoreIndices(subset_mult,&idxs_mult);
165: } else {
166: for (i=0;i<n;i++) {
167: gidxs[i] = leaf_data[i]-1;
168: }
169: }
170: ISCreateGeneral(PetscObjectComm((PetscObject)subset),n_n,gidxs,PETSC_OWN_POINTER,subset_n);
171: PetscFree2(leaf_data,root_data);
172: return(0);
173: }
176: /*@
177: ISCreateSubIS - Create a sub index set from a global index set selecting some components.
179: Collective on IS
181: Input Parmeters:
182: + is - the index set
183: - comps - which components we will extract from is
185: Output Parameters:
186: . subis - the new sub index set
188: Level: intermediate
190: Example usage:
191: We have an index set (is) living on 3 processes with the following values:
192: | 4 9 0 | 2 6 7 | 10 11 1|
193: and another index set (comps) used to indicate which components of is we want to take,
194: | 7 5 | 1 2 | 0 4|
195: The output index set (subis) should look like:
196: | 11 7 | 9 0 | 4 6|
198: .seealso: VecGetSubVector(), MatCreateSubMatrix()
199: @*/
200: PetscErrorCode ISCreateSubIS(IS is,IS comps,IS *subis)
201: {
202: PetscSF sf;
203: const PetscInt *is_indices,*comps_indices;
204: PetscInt *subis_indices,nroots,nleaves,*mine,i,lidx;
205: PetscMPIInt owner;
206: PetscSFNode *remote;
207: PetscErrorCode ierr;
208: MPI_Comm comm;
215: PetscObjectGetComm((PetscObject)is, &comm);
216: ISGetLocalSize(comps,&nleaves);
217: ISGetLocalSize(is,&nroots);
218: PetscMalloc1(nleaves,&remote);
219: PetscMalloc1(nleaves,&mine);
220: ISGetIndices(comps,&comps_indices);
221: /*
222: * Construct a PetscSF in which "is" data serves as roots and "subis" is leaves.
223: * Root data are sent to leaves using PetscSFBcast().
224: * */
225: for (i=0; i<nleaves; i++) {
226: mine[i] = i;
227: /* Connect a remote root with the current leaf. The value on the remote root
228: * will be received by the current local leaf.
229: * */
230: owner = -1;
231: lidx = -1;
232: PetscLayoutFindOwnerIndex(is->map,comps_indices[i],&owner,&lidx);
233: remote[i].rank = owner;
234: remote[i].index = lidx;
235: }
236: ISRestoreIndices(comps,&comps_indices);
237: PetscSFCreate(comm,&sf);
238: PetscSFSetFromOptions(sf);\
239: PetscSFSetGraph(sf,nroots,nleaves,mine,PETSC_OWN_POINTER,remote,PETSC_OWN_POINTER);
241: PetscMalloc1(nleaves,&subis_indices);
242: ISGetIndices(is, &is_indices);
243: PetscSFBcastBegin(sf,MPIU_INT,is_indices,subis_indices,MPI_REPLACE);
244: PetscSFBcastEnd(sf,MPIU_INT,is_indices,subis_indices,MPI_REPLACE);
245: ISRestoreIndices(is,&is_indices);
246: PetscSFDestroy(&sf);
247: ISCreateGeneral(comm,nleaves,subis_indices,PETSC_OWN_POINTER,subis);
248: return(0);
249: }
251: /*@
252: ISClearInfoCache - clear the cache of computed index set properties
254: Not collective
256: Input Parameters:
257: + is - the index set
258: - clear_permanent_local - whether to remove the permanent status of local properties
260: NOTE: because all processes must agree on the global permanent status of a property,
261: the permanent status can only be changed with ISSetInfo(), because this routine is not collective
263: Level: developer
265: .seealso: ISInfo, ISInfoType, ISSetInfo(), ISClearInfoCache()
267: @*/
268: PetscErrorCode ISClearInfoCache(IS is, PetscBool clear_permanent_local)
269: {
270: PetscInt i, j;
275: for (i = 0; i < IS_INFO_MAX; i++) {
276: if (clear_permanent_local) is->info_permanent[IS_LOCAL][i] = PETSC_FALSE;
277: for (j = 0; j < 2; j++) {
278: if (!is->info_permanent[j][i]) is->info[j][i] = IS_INFO_UNKNOWN;
279: }
280: }
281: return(0);
282: }
284: static PetscErrorCode ISSetInfo_Internal(IS is, ISInfo info, ISInfoType type, ISInfoBool ipermanent, PetscBool flg)
285: {
286: ISInfoBool iflg = flg ? IS_INFO_TRUE : IS_INFO_FALSE;
287: PetscInt itype = (type == IS_LOCAL) ? 0 : 1;
288: PetscBool permanent_set = (ipermanent == IS_INFO_UNKNOWN) ? PETSC_FALSE : PETSC_TRUE;
289: PetscBool permanent = (ipermanent == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
292: /* set this property */
293: is->info[itype][(int)info] = iflg;
294: if (permanent_set) is->info_permanent[itype][(int)info] = permanent;
295: /* set implications */
296: switch (info) {
297: case IS_SORTED:
298: if (flg && type == IS_GLOBAL) { /* an array that is globally sorted is also locally sorted */
299: is->info[IS_LOCAL][(int)info] = IS_INFO_TRUE;
300: /* global permanence implies local permanence */
301: if (permanent_set && permanent) is->info_permanent[IS_LOCAL][(int)info] = PETSC_TRUE;
302: }
303: if (!flg) { /* if an array is not sorted, it cannot be an interval or the identity */
304: is->info[itype][IS_INTERVAL] = IS_INFO_FALSE;
305: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
306: if (permanent_set) {
307: is->info_permanent[itype][IS_INTERVAL] = permanent;
308: is->info_permanent[itype][IS_IDENTITY] = permanent;
309: }
310: }
311: break;
312: case IS_UNIQUE:
313: if (flg && type == IS_GLOBAL) { /* an array that is globally unique is also locally unique */
314: is->info[IS_LOCAL][(int)info] = IS_INFO_TRUE;
315: /* global permanence implies local permanence */
316: if (permanent_set && permanent) is->info_permanent[IS_LOCAL][(int)info] = PETSC_TRUE;
317: }
318: if (!flg) { /* if an array is not unique, it cannot be a permutation, and interval, or the identity */
319: is->info[itype][IS_PERMUTATION] = IS_INFO_FALSE;
320: is->info[itype][IS_INTERVAL] = IS_INFO_FALSE;
321: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
322: if (permanent_set) {
323: is->info_permanent[itype][IS_PERMUTATION] = permanent;
324: is->info_permanent[itype][IS_INTERVAL] = permanent;
325: is->info_permanent[itype][IS_IDENTITY] = permanent;
326: }
327: }
328: break;
329: case IS_PERMUTATION:
330: if (flg) { /* an array that is a permutation is unique and is unique locally */
331: is->info[itype][IS_UNIQUE] = IS_INFO_TRUE;
332: is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
333: if (permanent_set && permanent) {
334: is->info_permanent[itype][IS_UNIQUE] = PETSC_TRUE;
335: is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
336: }
337: } else { /* an array that is not a permutation cannot be the identity */
338: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
339: if (permanent_set) is->info_permanent[itype][IS_IDENTITY] = permanent;
340: }
341: break;
342: case IS_INTERVAL:
343: if (flg) { /* an array that is an interval is sorted and unique */
344: is->info[itype][IS_SORTED] = IS_INFO_TRUE;
345: is->info[IS_LOCAL][IS_SORTED] = IS_INFO_TRUE;
346: is->info[itype][IS_UNIQUE] = IS_INFO_TRUE;
347: is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
348: if (permanent_set && permanent) {
349: is->info_permanent[itype][IS_SORTED] = PETSC_TRUE;
350: is->info_permanent[IS_LOCAL][IS_SORTED] = PETSC_TRUE;
351: is->info_permanent[itype][IS_UNIQUE] = PETSC_TRUE;
352: is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
353: }
354: } else { /* an array that is not an interval cannot be the identity */
355: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
356: if (permanent_set) is->info_permanent[itype][IS_IDENTITY] = permanent;
357: }
358: break;
359: case IS_IDENTITY:
360: if (flg) { /* an array that is the identity is sorted, unique, an interval, and a permutation */
361: is->info[itype][IS_SORTED] = IS_INFO_TRUE;
362: is->info[IS_LOCAL][IS_SORTED] = IS_INFO_TRUE;
363: is->info[itype][IS_UNIQUE] = IS_INFO_TRUE;
364: is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
365: is->info[itype][IS_PERMUTATION] = IS_INFO_TRUE;
366: is->info[itype][IS_INTERVAL] = IS_INFO_TRUE;
367: is->info[IS_LOCAL][IS_INTERVAL] = IS_INFO_TRUE;
368: if (permanent_set && permanent) {
369: is->info_permanent[itype][IS_SORTED] = PETSC_TRUE;
370: is->info_permanent[IS_LOCAL][IS_SORTED] = PETSC_TRUE;
371: is->info_permanent[itype][IS_UNIQUE] = PETSC_TRUE;
372: is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
373: is->info_permanent[itype][IS_PERMUTATION] = PETSC_TRUE;
374: is->info_permanent[itype][IS_INTERVAL] = PETSC_TRUE;
375: is->info_permanent[IS_LOCAL][IS_INTERVAL] = PETSC_TRUE;
376: }
377: }
378: break;
379: default:
380: if (type == IS_LOCAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
381: else SETERRQ(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
382: }
383: return(0);
384: }
386: /*@
387: ISSetInfo - Set known information about an index set.
389: Logically Collective on IS if type is IS_GLOBAL
391: Input Parameters:
392: + is - the index set
393: . info - describing a property of the index set, one of those listed below,
394: . type - IS_LOCAL if the information describes the local portion of the index set,
395: IS_GLOBAL if it describes the whole index set
396: . permanent - PETSC_TRUE if it is known that the property will persist through changes to the index set, PETSC_FALSE otherwise
397: If the user sets a property as permanently known, it will bypass computation of that property
398: - flg - set the described property as true (PETSC_TRUE) or false (PETSC_FALSE)
400: Info Describing IS Structure:
401: + IS_SORTED - the [local part of the] index set is sorted in ascending order
402: . IS_UNIQUE - each entry in the [local part of the] index set is unique
403: . IS_PERMUTATION - the [local part of the] index set is a permutation of the integers {0, 1, ..., N-1}, where N is the size of the [local part of the] index set
404: . IS_INTERVAL - the [local part of the] index set is equal to a contiguous range of integers {f, f + 1, ..., f + N-1}
405: - IS_IDENTITY - the [local part of the] index set is equal to the integers {0, 1, ..., N-1}
408: Notes:
409: If type is IS_GLOBAL, all processes that share the index set must pass the same value in flg
411: It is possible to set a property with ISSetInfo() that contradicts what would be previously computed with ISGetInfo()
413: Level: advanced
415: .seealso: ISInfo, ISInfoType, IS
417: @*/
418: PetscErrorCode ISSetInfo(IS is, ISInfo info, ISInfoType type, PetscBool permanent, PetscBool flg)
419: {
420: MPI_Comm comm, errcomm;
421: PetscMPIInt size;
427: comm = PetscObjectComm((PetscObject)is);
428: if (type == IS_GLOBAL) {
432: errcomm = comm;
433: } else {
434: errcomm = PETSC_COMM_SELF;
435: }
437: if (((int) info) <= IS_INFO_MIN || ((int) info) >= IS_INFO_MAX) SETERRQ1(errcomm,PETSC_ERR_ARG_OUTOFRANGE,"Options %d is out of range",(int)info);
439: MPI_Comm_size(comm, &size);
440: /* do not use global values if size == 1: it makes it easier to keep the implications straight */
441: if (size == 1) type = IS_LOCAL;
442: ISSetInfo_Internal(is, info, type, permanent ? IS_INFO_TRUE : IS_INFO_FALSE, flg);
443: return(0);
444: }
446: static PetscErrorCode ISGetInfo_Sorted(IS is, ISInfoType type, PetscBool *flg)
447: {
448: MPI_Comm comm;
449: PetscMPIInt size, rank;
453: comm = PetscObjectComm((PetscObject)is);
454: MPI_Comm_size(comm, &size);
455: MPI_Comm_size(comm, &rank);
456: if (type == IS_GLOBAL && is->ops->sortedglobal) {
457: (*is->ops->sortedglobal)(is,flg);
458: } else {
459: PetscBool sortedLocal = PETSC_FALSE;
461: /* determine if the array is locally sorted */
462: if (type == IS_GLOBAL && size > 1) {
463: /* call ISGetInfo so that a cached value will be used if possible */
464: ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sortedLocal);
465: } else if (is->ops->sortedlocal) {
466: (*is->ops->sortedlocal)(is,&sortedLocal);
467: } else {
468: /* default: get the local indices and directly check */
469: const PetscInt *idx;
470: PetscInt n;
472: ISGetIndices(is, &idx);
473: ISGetLocalSize(is, &n);
474: PetscSortedInt(n, idx, &sortedLocal);
475: ISRestoreIndices(is, &idx);
476: }
478: if (type == IS_LOCAL || size == 1) {
479: *flg = sortedLocal;
480: } else {
481: MPI_Allreduce(&sortedLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
482: if (*flg) {
483: PetscInt n, min = PETSC_MAX_INT, max = PETSC_MIN_INT;
484: PetscInt maxprev;
486: ISGetLocalSize(is, &n);
487: if (n) {ISGetMinMax(is, &min, &max);}
488: maxprev = PETSC_MIN_INT;
489: MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm);
490: if (rank && (maxprev > min)) sortedLocal = PETSC_FALSE;
491: MPI_Allreduce(&sortedLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
492: }
493: }
494: }
495: return(0);
496: }
498: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[]);
500: static PetscErrorCode ISGetInfo_Unique(IS is, ISInfoType type, PetscBool *flg)
501: {
502: MPI_Comm comm;
503: PetscMPIInt size, rank;
504: PetscInt i;
508: comm = PetscObjectComm((PetscObject)is);
509: MPI_Comm_size(comm, &size);
510: MPI_Comm_size(comm, &rank);
511: if (type == IS_GLOBAL && is->ops->uniqueglobal) {
512: (*is->ops->uniqueglobal)(is,flg);
513: } else {
514: PetscBool uniqueLocal;
515: PetscInt n = -1;
516: PetscInt *idx = NULL;
518: /* determine if the array is locally unique */
519: if (type == IS_GLOBAL && size > 1) {
520: /* call ISGetInfo so that a cached value will be used if possible */
521: ISGetInfo(is, IS_UNIQUE, IS_LOCAL, PETSC_TRUE, &uniqueLocal);
522: } else if (is->ops->uniquelocal) {
523: (*is->ops->uniquelocal)(is,&uniqueLocal);
524: } else {
525: /* default: get the local indices and directly check */
526: uniqueLocal = PETSC_TRUE;
527: ISGetLocalSize(is, &n);
528: PetscMalloc1(n, &idx);
529: ISGetIndicesCopy(is, idx);
530: PetscIntSortSemiOrdered(n, idx);
531: for (i = 1; i < n; i++) if (idx[i] == idx[i-1]) break;
532: if (i < n) uniqueLocal = PETSC_FALSE;
533: }
535: PetscFree(idx);
536: if (type == IS_LOCAL || size == 1) {
537: *flg = uniqueLocal;
538: } else {
539: MPI_Allreduce(&uniqueLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
540: if (*flg) {
541: PetscInt min = PETSC_MAX_INT, max = PETSC_MIN_INT, maxprev;
543: if (!idx) {
544: ISGetLocalSize(is, &n);
545: PetscMalloc1(n, &idx);
546: ISGetIndicesCopy(is, idx);
547: }
548: PetscParallelSortInt(is->map, is->map, idx, idx);
549: if (n) {
550: min = idx[0];
551: max = idx[n - 1];
552: }
553: for (i = 1; i < n; i++) if (idx[i] == idx[i-1]) break;
554: if (i < n) uniqueLocal = PETSC_FALSE;
555: maxprev = PETSC_MIN_INT;
556: MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm);
557: if (rank && (maxprev == min)) uniqueLocal = PETSC_FALSE;
558: MPI_Allreduce(&uniqueLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
559: }
560: }
561: PetscFree(idx);
562: }
563: return(0);
564: }
566: static PetscErrorCode ISGetInfo_Permutation(IS is, ISInfoType type, PetscBool *flg)
567: {
568: MPI_Comm comm;
569: PetscMPIInt size, rank;
573: comm = PetscObjectComm((PetscObject)is);
574: MPI_Comm_size(comm, &size);
575: MPI_Comm_size(comm, &rank);
576: if (type == IS_GLOBAL && is->ops->permglobal) {
577: (*is->ops->permglobal)(is,flg);
578: } else if (type == IS_LOCAL && is->ops->permlocal) {
579: (*is->ops->permlocal)(is,flg);
580: } else {
581: PetscBool permLocal;
582: PetscInt n, i, rStart;
583: PetscInt *idx;
585: ISGetLocalSize(is, &n);
586: PetscMalloc1(n, &idx);
587: ISGetIndicesCopy(is, idx);
588: if (type == IS_GLOBAL) {
589: PetscParallelSortInt(is->map, is->map, idx, idx);
590: PetscLayoutGetRange(is->map, &rStart, NULL);
591: } else {
592: PetscIntSortSemiOrdered(n, idx);
593: rStart = 0;
594: }
595: permLocal = PETSC_TRUE;
596: for (i = 0; i < n; i++) {
597: if (idx[i] != rStart + i) break;
598: }
599: if (i < n) permLocal = PETSC_FALSE;
600: if (type == IS_LOCAL || size == 1) {
601: *flg = permLocal;
602: } else {
603: MPI_Allreduce(&permLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
604: }
605: PetscFree(idx);
606: }
607: return(0);
608: }
610: static PetscErrorCode ISGetInfo_Interval(IS is, ISInfoType type, PetscBool *flg)
611: {
612: MPI_Comm comm;
613: PetscMPIInt size, rank;
614: PetscInt i;
618: comm = PetscObjectComm((PetscObject)is);
619: MPI_Comm_size(comm, &size);
620: MPI_Comm_size(comm, &rank);
621: if (type == IS_GLOBAL && is->ops->intervalglobal) {
622: (*is->ops->intervalglobal)(is,flg);
623: } else {
624: PetscBool intervalLocal;
626: /* determine if the array is locally an interval */
627: if (type == IS_GLOBAL && size > 1) {
628: /* call ISGetInfo so that a cached value will be used if possible */
629: ISGetInfo(is, IS_INTERVAL, IS_LOCAL, PETSC_TRUE, &intervalLocal);
630: } else if (is->ops->intervallocal) {
631: (*is->ops->intervallocal)(is,&intervalLocal);
632: } else {
633: PetscInt n;
634: const PetscInt *idx;
635: /* default: get the local indices and directly check */
636: intervalLocal = PETSC_TRUE;
637: ISGetLocalSize(is, &n);
638: PetscMalloc1(n, &idx);
639: ISGetIndices(is, &idx);
640: for (i = 1; i < n; i++) if (idx[i] != idx[i-1] + 1) break;
641: if (i < n) intervalLocal = PETSC_FALSE;
642: ISRestoreIndices(is, &idx);
643: }
645: if (type == IS_LOCAL || size == 1) {
646: *flg = intervalLocal;
647: } else {
648: MPI_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
649: if (*flg) {
650: PetscInt n, min = PETSC_MAX_INT, max = PETSC_MIN_INT;
651: PetscInt maxprev;
653: ISGetLocalSize(is, &n);
654: if (n) {ISGetMinMax(is, &min, &max);}
655: maxprev = PETSC_MIN_INT;
656: MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm);
657: if (rank && n && (maxprev != min - 1)) intervalLocal = PETSC_FALSE;
658: MPI_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
659: }
660: }
661: }
662: return(0);
663: }
665: static PetscErrorCode ISGetInfo_Identity(IS is, ISInfoType type, PetscBool *flg)
666: {
667: MPI_Comm comm;
668: PetscMPIInt size, rank;
672: comm = PetscObjectComm((PetscObject)is);
673: MPI_Comm_size(comm, &size);
674: MPI_Comm_size(comm, &rank);
675: if (type == IS_GLOBAL && is->ops->intervalglobal) {
676: PetscBool isinterval;
678: (*is->ops->intervalglobal)(is,&isinterval);
679: *flg = PETSC_FALSE;
680: if (isinterval) {
681: PetscInt min;
683: ISGetMinMax(is, &min, NULL);
684: MPI_Bcast(&min, 1, MPIU_INT, 0, comm);
685: if (min == 0) *flg = PETSC_TRUE;
686: }
687: } else if (type == IS_LOCAL && is->ops->intervallocal) {
688: PetscBool isinterval;
690: (*is->ops->intervallocal)(is,&isinterval);
691: *flg = PETSC_FALSE;
692: if (isinterval) {
693: PetscInt min;
695: ISGetMinMax(is, &min, NULL);
696: if (min == 0) *flg = PETSC_TRUE;
697: }
698: } else {
699: PetscBool identLocal;
700: PetscInt n, i, rStart;
701: const PetscInt *idx;
703: ISGetLocalSize(is, &n);
704: ISGetIndices(is, &idx);
705: PetscLayoutGetRange(is->map, &rStart, NULL);
706: identLocal = PETSC_TRUE;
707: for (i = 0; i < n; i++) {
708: if (idx[i] != rStart + i) break;
709: }
710: if (i < n) identLocal = PETSC_FALSE;
711: if (type == IS_LOCAL || size == 1) {
712: *flg = identLocal;
713: } else {
714: MPI_Allreduce(&identLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
715: }
716: ISRestoreIndices(is, &idx);
717: }
718: return(0);
719: }
721: /*@
722: ISGetInfo - Determine whether an index set satisfies a given property
724: Collective or logically collective on IS if the type is IS_GLOBAL (logically collective if the value of the property has been permanently set with ISSetInfo())
726: Input Parameters:
727: + is - the index set
728: . info - describing a property of the index set, one of those listed in the documentation of ISSetInfo()
729: . compute - if PETSC_FALSE, the property will not be computed if it is not already known and the property will be assumed to be false
730: - type - whether the property is local (IS_LOCAL) or global (IS_GLOBAL)
732: Output Parameter:
733: . flg - wheter the property is true (PETSC_TRUE) or false (PETSC_FALSE)
735: Note: ISGetInfo uses cached values when possible, which will be incorrect if ISSetInfo() has been called with incorrect information. To clear cached values, use ISClearInfoCache().
737: Level: advanced
739: .seealso: ISInfo, ISInfoType, ISSetInfo(), ISClearInfoCache()
741: @*/
742: PetscErrorCode ISGetInfo(IS is, ISInfo info, ISInfoType type, PetscBool compute, PetscBool *flg)
743: {
744: MPI_Comm comm, errcomm;
745: PetscMPIInt rank, size;
746: PetscInt itype;
747: PetscBool hasprop;
748: PetscBool infer;
754: comm = PetscObjectComm((PetscObject)is);
755: if (type == IS_GLOBAL) {
757: errcomm = comm;
758: } else {
759: errcomm = PETSC_COMM_SELF;
760: }
762: MPI_Comm_size(comm, &size);
763: MPI_Comm_rank(comm, &rank);
765: if (((int) info) <= IS_INFO_MIN || ((int) info) >= IS_INFO_MAX) SETERRQ1(errcomm,PETSC_ERR_ARG_OUTOFRANGE,"Options %d is out of range",(int)info);
766: if (size == 1) type = IS_LOCAL;
767: itype = (type == IS_LOCAL) ? 0 : 1;
768: hasprop = PETSC_FALSE;
769: infer = PETSC_FALSE;
770: if (is->info_permanent[itype][(int)info]) {
771: hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
772: infer = PETSC_TRUE;
773: } else if ((itype == IS_LOCAL) && (is->info[IS_LOCAL][info] != IS_INFO_UNKNOWN)) {
774: /* we can cache local properties as long as we clear them when the IS changes */
775: /* NOTE: we only cache local values because there is no ISAssemblyBegin()/ISAssemblyEnd(),
776: so we have no way of knowing when a cached value has been invalidated by changes on a different process */
777: hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
778: infer = PETSC_TRUE;
779: } else if (compute) {
780: switch (info) {
781: case IS_SORTED:
782: ISGetInfo_Sorted(is, type, &hasprop);
783: break;
784: case IS_UNIQUE:
785: ISGetInfo_Unique(is, type, &hasprop);
786: break;
787: case IS_PERMUTATION:
788: ISGetInfo_Permutation(is, type, &hasprop);
789: break;
790: case IS_INTERVAL:
791: ISGetInfo_Interval(is, type, &hasprop);
792: break;
793: case IS_IDENTITY:
794: ISGetInfo_Identity(is, type, &hasprop);
795: break;
796: default:
797: SETERRQ(errcomm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
798: }
799: infer = PETSC_TRUE;
800: }
801: /* call ISSetInfo_Internal to keep all of the implications straight */
802: if (infer) {ISSetInfo_Internal(is, info, type, IS_INFO_UNKNOWN, hasprop);}
803: *flg = hasprop;
804: return(0);
805: }
807: static PetscErrorCode ISCopyInfo(IS source, IS dest)
808: {
812: PetscArraycpy(&dest->info[0], &source->info[0], 2);
813: PetscArraycpy(&dest->info_permanent[0], &source->info_permanent[0], 2);
814: return(0);
815: }
817: /*@
818: ISIdentity - Determines whether index set is the identity mapping.
820: Collective on IS
822: Input Parmeters:
823: . is - the index set
825: Output Parameters:
826: . ident - PETSC_TRUE if an identity, else PETSC_FALSE
828: Level: intermediate
830: Note: If ISSetIdentity() (or ISSetInfo() for a permanent property) has been called,
831: ISIdentity() will return its answer without communication between processes, but
832: otherwise the output ident will be computed from ISGetInfo(),
833: which may require synchronization on the communicator of IS. To avoid this computation,
834: call ISGetInfo() directly with the compute flag set to PETSC_FALSE, and ident will be assumed false.
836: .seealso: ISSetIdentity(), ISGetInfo()
837: @*/
838: PetscErrorCode ISIdentity(IS is,PetscBool *ident)
839: {
845: ISGetInfo(is,IS_IDENTITY,IS_GLOBAL,PETSC_TRUE,ident);
846: return(0);
847: }
849: /*@
850: ISSetIdentity - Informs the index set that it is an identity.
852: Logically Collective on IS
854: Input Parmeters:
855: . is - the index set
857: Level: intermediate
859: Note: The IS will be considered the identity permanently, even if indices have been changes (for example, with
860: ISGeneralSetIndices()). It's a good idea to only set this property if the IS will not change in the future.
861: To clear this property, use ISClearInfoCache().
863: .seealso: ISIdentity(), ISSetInfo(), ISClearInfoCache()
864: @*/
865: PetscErrorCode ISSetIdentity(IS is)
866: {
871: ISSetInfo(is,IS_IDENTITY,IS_GLOBAL,PETSC_TRUE,PETSC_TRUE);
872: return(0);
873: }
875: /*@
876: ISContiguousLocal - Locates an index set with contiguous range within a global range, if possible
878: Not Collective
880: Input Parmeters:
881: + is - the index set
882: . gstart - global start
883: - gend - global end
885: Output Parameters:
886: + start - start of contiguous block, as an offset from gstart
887: - contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
889: Level: developer
891: .seealso: ISGetLocalSize(), VecGetOwnershipRange()
892: @*/
893: PetscErrorCode ISContiguousLocal(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
894: {
901: *start = -1;
902: *contig = PETSC_FALSE;
903: if (is->ops->contiguous) {
904: (*is->ops->contiguous)(is,gstart,gend,start,contig);
905: }
906: return(0);
907: }
909: /*@
910: ISPermutation - PETSC_TRUE or PETSC_FALSE depending on whether the
911: index set has been declared to be a permutation.
913: Logically Collective on IS
915: Input Parmeters:
916: . is - the index set
918: Output Parameters:
919: . perm - PETSC_TRUE if a permutation, else PETSC_FALSE
921: Level: intermediate
923: Note: If it is not alread known that the IS is a permutation (if ISSetPermutation()
924: or ISSetInfo() has not been called), this routine will not attempt to compute
925: whether the index set is a permutation and will assume perm is PETSC_FALSE.
926: To compute the value when it is not already known, use ISGetInfo() with
927: the compute flag set to PETSC_TRUE.
929: .seealso: ISSetPermutation(), ISGetInfo()
930: @*/
931: PetscErrorCode ISPermutation(IS is,PetscBool *perm)
932: {
938: ISGetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_FALSE,perm);
939: return(0);
940: }
942: /*@
943: ISSetPermutation - Informs the index set that it is a permutation.
945: Logically Collective on IS
947: Input Parmeters:
948: . is - the index set
950: Level: intermediate
953: The debug version of the libraries (./configure --with-debugging=1) checks if the
954: index set is actually a permutation. The optimized version just believes you.
956: Note: The IS will be considered a permutation permanently, even if indices have been changes (for example, with
957: ISGeneralSetIndices()). It's a good idea to only set this property if the IS will not change in the future.
958: To clear this property, use ISClearInfoCache().
960: .seealso: ISPermutation(), ISSetInfo(), ISClearInfoCache().
961: @*/
962: PetscErrorCode ISSetPermutation(IS is)
963: {
968: if (PetscDefined(USE_DEBUG)) {
969: PetscMPIInt size;
971: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
972: if (size == 1) {
973: PetscInt i,n,*idx;
974: const PetscInt *iidx;
976: ISGetSize(is,&n);
977: PetscMalloc1(n,&idx);
978: ISGetIndices(is,&iidx);
979: PetscArraycpy(idx,iidx,n);
980: PetscIntSortSemiOrdered(n,idx);
981: for (i=0; i<n; i++) {
982: if (idx[i] != i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Index set is not a permutation");
983: }
984: PetscFree(idx);
985: ISRestoreIndices(is,&iidx);
986: }
987: }
988: ISSetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_TRUE,PETSC_TRUE);
989: return(0);
990: }
992: /*@C
993: ISDestroy - Destroys an index set.
995: Collective on IS
997: Input Parameters:
998: . is - the index set
1000: Level: beginner
1002: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlocked()
1003: @*/
1004: PetscErrorCode ISDestroy(IS *is)
1005: {
1009: if (!*is) return(0);
1011: if (--((PetscObject)(*is))->refct > 0) {*is = NULL; return(0);}
1012: if ((*is)->complement) {
1013: PetscInt refcnt;
1014: PetscObjectGetReference((PetscObject)((*is)->complement), &refcnt);
1015: if (refcnt > 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Nonlocal IS has not been restored");
1016: ISDestroy(&(*is)->complement);
1017: }
1018: if ((*is)->ops->destroy) {
1019: (*(*is)->ops->destroy)(*is);
1020: }
1021: PetscLayoutDestroy(&(*is)->map);
1022: /* Destroy local representations of offproc data. */
1023: PetscFree((*is)->total);
1024: PetscFree((*is)->nonlocal);
1025: PetscHeaderDestroy(is);
1026: return(0);
1027: }
1029: /*@
1030: ISInvertPermutation - Creates a new permutation that is the inverse of
1031: a given permutation.
1033: Collective on IS
1035: Input Parameter:
1036: + is - the index set
1037: - nlocal - number of indices on this processor in result (ignored for 1 proccessor) or
1038: use PETSC_DECIDE
1040: Output Parameter:
1041: . isout - the inverse permutation
1043: Level: intermediate
1045: Notes:
1046: For parallel index sets this does the complete parallel permutation, but the
1047: code is not efficient for huge index sets (10,000,000 indices).
1049: @*/
1050: PetscErrorCode ISInvertPermutation(IS is,PetscInt nlocal,IS *isout)
1051: {
1052: PetscBool isperm, isidentity, issame;
1058: ISGetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_TRUE,&isperm);
1059: if (!isperm) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_WRONG,"Not a permutation");
1060: ISGetInfo(is,IS_IDENTITY,IS_GLOBAL,PETSC_TRUE,&isidentity);
1061: issame = PETSC_FALSE;
1062: if (isidentity) {
1063: PetscInt n;
1064: PetscBool isallsame;
1066: ISGetLocalSize(is, &n);
1067: issame = (PetscBool) (n == nlocal);
1068: MPI_Allreduce(&issame, &isallsame, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is));
1069: issame = isallsame;
1070: }
1071: if (issame) {
1072: ISDuplicate(is,isout);
1073: } else {
1074: (*is->ops->invertpermutation)(is,nlocal,isout);
1075: ISSetPermutation(*isout);
1076: }
1077: return(0);
1078: }
1080: /*@
1081: ISGetSize - Returns the global length of an index set.
1083: Not Collective
1085: Input Parameter:
1086: . is - the index set
1088: Output Parameter:
1089: . size - the global size
1091: Level: beginner
1094: @*/
1095: PetscErrorCode ISGetSize(IS is,PetscInt *size)
1096: {
1100: *size = is->map->N;
1101: return(0);
1102: }
1104: /*@
1105: ISGetLocalSize - Returns the local (processor) length of an index set.
1107: Not Collective
1109: Input Parameter:
1110: . is - the index set
1112: Output Parameter:
1113: . size - the local size
1115: Level: beginner
1117: @*/
1118: PetscErrorCode ISGetLocalSize(IS is,PetscInt *size)
1119: {
1123: *size = is->map->n;
1124: return(0);
1125: }
1127: /*@
1128: ISGetLayout - get PetscLayout describing index set layout
1130: Not Collective
1132: Input Arguments:
1133: . is - the index set
1135: Output Arguments:
1136: . map - the layout
1138: Level: developer
1140: .seealso: ISGetSize(), ISGetLocalSize()
1141: @*/
1142: PetscErrorCode ISGetLayout(IS is,PetscLayout *map)
1143: {
1148: *map = is->map;
1149: return(0);
1150: }
1152: /*@C
1153: ISGetIndices - Returns a pointer to the indices. The user should call
1154: ISRestoreIndices() after having looked at the indices. The user should
1155: NOT change the indices.
1157: Not Collective
1159: Input Parameter:
1160: . is - the index set
1162: Output Parameter:
1163: . ptr - the location to put the pointer to the indices
1165: Fortran Note:
1166: This routine has two different interfaces from Fortran; the first is not recommend, it does not require Fortran 90
1167: $ IS is
1168: $ integer is_array(1)
1169: $ PetscOffset i_is
1170: $ int ierr
1171: $ call ISGetIndices(is,is_array,i_is,ierr)
1172: $
1173: $ Access first local entry in list
1174: $ value = is_array(i_is + 1)
1175: $
1176: $ ...... other code
1177: $ call ISRestoreIndices(is,is_array,i_is,ierr)
1178: The second Fortran interface is recommended.
1179: $ use petscisdef
1180: $ PetscInt, pointer :: array(:)
1181: $ PetscErrorCode ierr
1182: $ IS i
1183: $ call ISGetIndicesF90(i,array,ierr)
1187: See the Fortran chapter of the users manual and
1188: petsc/src/is/[tutorials,tests] for details.
1190: Level: intermediate
1193: .seealso: ISRestoreIndices(), ISGetIndicesF90()
1194: @*/
1195: PetscErrorCode ISGetIndices(IS is,const PetscInt *ptr[])
1196: {
1202: (*is->ops->getindices)(is,ptr);
1203: return(0);
1204: }
1206: /*@C
1207: ISGetMinMax - Gets the minimum and maximum values in an IS
1209: Not Collective
1211: Input Parameter:
1212: . is - the index set
1214: Output Parameter:
1215: + min - the minimum value
1216: - max - the maximum value
1218: Level: intermediate
1220: Notes:
1221: Empty index sets return min=PETSC_MAX_INT and max=PETSC_MIN_INT.
1222: In parallel, it returns the min and max of the local portion of the IS
1225: .seealso: ISGetIndices(), ISRestoreIndices(), ISGetIndicesF90()
1226: @*/
1227: PetscErrorCode ISGetMinMax(IS is,PetscInt *min,PetscInt *max)
1228: {
1231: if (min) *min = is->min;
1232: if (max) *max = is->max;
1233: return(0);
1234: }
1236: /*@
1237: ISLocate - determine the location of an index within the local component of an index set
1239: Not Collective
1241: Input Parameter:
1242: + is - the index set
1243: - key - the search key
1245: Output Parameter:
1246: . location - if >= 0, a location within the index set that is equal to the key, otherwise the key is not in the index set
1248: Level: intermediate
1249: @*/
1250: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
1251: {
1255: if (is->ops->locate) {
1256: (*is->ops->locate)(is,key,location);
1257: } else {
1258: PetscInt numIdx;
1259: PetscBool sorted;
1260: const PetscInt *idx;
1262: ISGetLocalSize(is,&numIdx);
1263: ISGetIndices(is,&idx);
1264: ISSorted(is,&sorted);
1265: if (sorted) {
1266: PetscFindInt(key,numIdx,idx,location);
1267: } else {
1268: PetscInt i;
1270: *location = -1;
1271: for (i = 0; i < numIdx; i++) {
1272: if (idx[i] == key) {
1273: *location = i;
1274: break;
1275: }
1276: }
1277: }
1278: ISRestoreIndices(is,&idx);
1279: }
1280: return(0);
1281: }
1283: /*@C
1284: ISRestoreIndices - Restores an index set to a usable state after a call
1285: to ISGetIndices().
1287: Not Collective
1289: Input Parameters:
1290: + is - the index set
1291: - ptr - the pointer obtained by ISGetIndices()
1293: Fortran Note:
1294: This routine is used differently from Fortran
1295: $ IS is
1296: $ integer is_array(1)
1297: $ PetscOffset i_is
1298: $ int ierr
1299: $ call ISGetIndices(is,is_array,i_is,ierr)
1300: $
1301: $ Access first local entry in list
1302: $ value = is_array(i_is + 1)
1303: $
1304: $ ...... other code
1305: $ call ISRestoreIndices(is,is_array,i_is,ierr)
1307: See the Fortran chapter of the users manual and
1308: petsc/src/vec/is/tests for details.
1310: Level: intermediate
1312: Note:
1313: This routine zeros out ptr. This is to prevent accidental us of the array after it has been restored.
1315: .seealso: ISGetIndices(), ISRestoreIndicesF90()
1316: @*/
1317: PetscErrorCode ISRestoreIndices(IS is,const PetscInt *ptr[])
1318: {
1324: if (is->ops->restoreindices) {
1325: (*is->ops->restoreindices)(is,ptr);
1326: }
1327: return(0);
1328: }
1330: static PetscErrorCode ISGatherTotal_Private(IS is)
1331: {
1333: PetscInt i,n,N;
1334: const PetscInt *lindices;
1335: MPI_Comm comm;
1336: PetscMPIInt rank,size,*sizes = NULL,*offsets = NULL,nn;
1341: PetscObjectGetComm((PetscObject)is,&comm);
1342: MPI_Comm_size(comm,&size);
1343: MPI_Comm_rank(comm,&rank);
1344: ISGetLocalSize(is,&n);
1345: PetscMalloc2(size,&sizes,size,&offsets);
1347: PetscMPIIntCast(n,&nn);
1348: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
1349: offsets[0] = 0;
1350: for (i=1; i<size; ++i) offsets[i] = offsets[i-1] + sizes[i-1];
1351: N = offsets[size-1] + sizes[size-1];
1353: PetscMalloc1(N,&(is->total));
1354: ISGetIndices(is,&lindices);
1355: MPI_Allgatherv((void*)lindices,nn,MPIU_INT,is->total,sizes,offsets,MPIU_INT,comm);
1356: ISRestoreIndices(is,&lindices);
1357: is->local_offset = offsets[rank];
1358: PetscFree2(sizes,offsets);
1359: return(0);
1360: }
1362: /*@C
1363: ISGetTotalIndices - Retrieve an array containing all indices across the communicator.
1365: Collective on IS
1367: Input Parameter:
1368: . is - the index set
1370: Output Parameter:
1371: . indices - total indices with rank 0 indices first, and so on; total array size is
1372: the same as returned with ISGetSize().
1374: Level: intermediate
1376: Notes:
1377: this is potentially nonscalable, but depends on the size of the total index set
1378: and the size of the communicator. This may be feasible for index sets defined on
1379: subcommunicators, such that the set size does not grow with PETSC_WORLD_COMM.
1380: Note also that there is no way to tell where the local part of the indices starts
1381: (use ISGetIndices() and ISGetNonlocalIndices() to retrieve just the local and just
1382: the nonlocal part (complement), respectively).
1384: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices(), ISGetSize()
1385: @*/
1386: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
1387: {
1389: PetscMPIInt size;
1394: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
1395: if (size == 1) {
1396: (*is->ops->getindices)(is,indices);
1397: } else {
1398: if (!is->total) {
1399: ISGatherTotal_Private(is);
1400: }
1401: *indices = is->total;
1402: }
1403: return(0);
1404: }
1406: /*@C
1407: ISRestoreTotalIndices - Restore the index array obtained with ISGetTotalIndices().
1409: Not Collective.
1411: Input Parameter:
1412: + is - the index set
1413: - indices - index array; must be the array obtained with ISGetTotalIndices()
1415: Level: intermediate
1417: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices()
1418: @*/
1419: PetscErrorCode ISRestoreTotalIndices(IS is, const PetscInt *indices[])
1420: {
1422: PetscMPIInt size;
1427: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
1428: if (size == 1) {
1429: (*is->ops->restoreindices)(is,indices);
1430: } else {
1431: 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.");
1432: }
1433: return(0);
1434: }
1435: /*@C
1436: ISGetNonlocalIndices - Retrieve an array of indices from remote processors
1437: in this communicator.
1439: Collective on IS
1441: Input Parameter:
1442: . is - the index set
1444: Output Parameter:
1445: . indices - indices with rank 0 indices first, and so on, omitting
1446: the current rank. Total number of indices is the difference
1447: total and local, obtained with ISGetSize() and ISGetLocalSize(),
1448: respectively.
1450: Level: intermediate
1452: Notes:
1453: restore the indices using ISRestoreNonlocalIndices().
1454: The same scalability considerations as those for ISGetTotalIndices
1455: apply here.
1457: .seealso: ISGetTotalIndices(), ISRestoreNonlocalIndices(), ISGetSize(), ISGetLocalSize().
1458: @*/
1459: PetscErrorCode ISGetNonlocalIndices(IS is, const PetscInt *indices[])
1460: {
1462: PetscMPIInt size;
1463: PetscInt n, N;
1468: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
1469: if (size == 1) *indices = NULL;
1470: else {
1471: if (!is->total) {
1472: ISGatherTotal_Private(is);
1473: }
1474: ISGetLocalSize(is,&n);
1475: ISGetSize(is,&N);
1476: PetscMalloc1(N-n, &(is->nonlocal));
1477: PetscArraycpy(is->nonlocal, is->total, is->local_offset);
1478: PetscArraycpy(is->nonlocal+is->local_offset, is->total+is->local_offset+n,N - is->local_offset - n);
1479: *indices = is->nonlocal;
1480: }
1481: return(0);
1482: }
1484: /*@C
1485: ISRestoreTotalIndices - Restore the index array obtained with ISGetNonlocalIndices().
1487: Not Collective.
1489: Input Parameter:
1490: + is - the index set
1491: - indices - index array; must be the array obtained with ISGetNonlocalIndices()
1493: Level: intermediate
1495: .seealso: ISGetTotalIndices(), ISGetNonlocalIndices(), ISRestoreTotalIndices()
1496: @*/
1497: PetscErrorCode ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
1498: {
1502: 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.");
1503: return(0);
1504: }
1506: /*@
1507: ISGetNonlocalIS - Gather all nonlocal indices for this IS and present
1508: them as another sequential index set.
1511: Collective on IS
1513: Input Parameter:
1514: . is - the index set
1516: Output Parameter:
1517: . complement - sequential IS with indices identical to the result of
1518: ISGetNonlocalIndices()
1520: Level: intermediate
1522: Notes:
1523: complement represents the result of ISGetNonlocalIndices as an IS.
1524: Therefore scalability issues similar to ISGetNonlocalIndices apply.
1525: The resulting IS must be restored using ISRestoreNonlocalIS().
1527: .seealso: ISGetNonlocalIndices(), ISRestoreNonlocalIndices(), ISAllGather(), ISGetSize()
1528: @*/
1529: PetscErrorCode ISGetNonlocalIS(IS is, IS *complement)
1530: {
1536: /* Check if the complement exists already. */
1537: if (is->complement) {
1538: *complement = is->complement;
1539: PetscObjectReference((PetscObject)(is->complement));
1540: } else {
1541: PetscInt N, n;
1542: const PetscInt *idx;
1543: ISGetSize(is, &N);
1544: ISGetLocalSize(is,&n);
1545: ISGetNonlocalIndices(is, &idx);
1546: ISCreateGeneral(PETSC_COMM_SELF, N-n,idx, PETSC_USE_POINTER, &(is->complement));
1547: PetscObjectReference((PetscObject)is->complement);
1548: *complement = is->complement;
1549: }
1550: return(0);
1551: }
1554: /*@
1555: ISRestoreNonlocalIS - Restore the IS obtained with ISGetNonlocalIS().
1557: Not collective.
1559: Input Parameter:
1560: + is - the index set
1561: - complement - index set of is's nonlocal indices
1563: Level: intermediate
1566: .seealso: ISGetNonlocalIS(), ISGetNonlocalIndices(), ISRestoreNonlocalIndices()
1567: @*/
1568: PetscErrorCode ISRestoreNonlocalIS(IS is, IS *complement)
1569: {
1571: PetscInt refcnt;
1576: if (*complement != is->complement) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Complement IS being restored was not obtained with ISGetNonlocalIS()");
1577: PetscObjectGetReference((PetscObject)(is->complement), &refcnt);
1578: if (refcnt <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Duplicate call to ISRestoreNonlocalIS() detected");
1579: PetscObjectDereference((PetscObject)(is->complement));
1580: return(0);
1581: }
1583: /*@C
1584: ISViewFromOptions - View from Options
1586: Collective on IS
1588: Input Parameters:
1589: + A - the index set
1590: . obj - Optional object
1591: - name - command line option
1593: Level: intermediate
1594: .seealso: IS, ISView, PetscObjectViewFromOptions(), ISCreate()
1595: @*/
1596: PetscErrorCode ISViewFromOptions(IS A,PetscObject obj,const char name[])
1597: {
1602: PetscObjectViewFromOptions((PetscObject)A,obj,name);
1603: return(0);
1604: }
1606: /*@C
1607: ISView - Displays an index set.
1609: Collective on IS
1611: Input Parameters:
1612: + is - the index set
1613: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
1615: Level: intermediate
1617: .seealso: PetscViewerASCIIOpen()
1618: @*/
1619: PetscErrorCode ISView(IS is,PetscViewer viewer)
1620: {
1625: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is),&viewer);}
1629: PetscObjectPrintClassNamePrefixType((PetscObject)is,viewer);
1630: PetscLogEventBegin(IS_View,is,viewer,0,0);
1631: (*is->ops->view)(is,viewer);
1632: PetscLogEventEnd(IS_View,is,viewer,0,0);
1633: return(0);
1634: }
1636: /*@
1637: ISLoad - Loads a vector that has been stored in binary or HDF5 format with ISView().
1639: Collective on PetscViewer
1641: Input Parameters:
1642: + is - the newly loaded vector, this needs to have been created with ISCreate() or some related function before a call to ISLoad().
1643: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or HDF5 file viewer, obtained from PetscViewerHDF5Open()
1645: Level: intermediate
1647: Notes:
1648: IF using HDF5, you must assign the IS the same name as was used in the IS
1649: that was stored in the file using PetscObjectSetName(). Otherwise you will
1650: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"
1652: .seealso: PetscViewerBinaryOpen(), ISView(), MatLoad(), VecLoad()
1653: @*/
1654: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
1655: {
1656: PetscBool isbinary, ishdf5;
1663: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
1664: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
1665: if (!isbinary && !ishdf5) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1666: if (!((PetscObject)is)->type_name) {ISSetType(is, ISGENERAL);}
1667: PetscLogEventBegin(IS_Load,is,viewer,0,0);
1668: (*is->ops->load)(is, viewer);
1669: PetscLogEventEnd(IS_Load,is,viewer,0,0);
1670: return(0);
1671: }
1673: /*@
1674: ISSort - Sorts the indices of an index set.
1676: Collective on IS
1678: Input Parameters:
1679: . is - the index set
1681: Level: intermediate
1684: .seealso: ISSortRemoveDups(), ISSorted()
1685: @*/
1686: PetscErrorCode ISSort(IS is)
1687: {
1692: (*is->ops->sort)(is);
1693: ISSetInfo(is,IS_SORTED,IS_LOCAL,is->info_permanent[IS_LOCAL][IS_SORTED],PETSC_TRUE);
1694: return(0);
1695: }
1697: /*@
1698: ISSortRemoveDups - Sorts the indices of an index set, removing duplicates.
1700: Collective on IS
1702: Input Parameters:
1703: . is - the index set
1705: Level: intermediate
1708: .seealso: ISSort(), ISSorted()
1709: @*/
1710: PetscErrorCode ISSortRemoveDups(IS is)
1711: {
1716: ISClearInfoCache(is,PETSC_FALSE);
1717: (*is->ops->sortremovedups)(is);
1718: ISSetInfo(is,IS_SORTED,IS_LOCAL,is->info_permanent[IS_LOCAL][IS_SORTED],PETSC_TRUE);
1719: ISSetInfo(is,IS_UNIQUE,IS_LOCAL,is->info_permanent[IS_LOCAL][IS_UNIQUE],PETSC_TRUE);
1720: return(0);
1721: }
1723: /*@
1724: ISToGeneral - Converts an IS object of any type to ISGENERAL type
1726: Collective on IS
1728: Input Parameters:
1729: . is - the index set
1731: Level: intermediate
1734: .seealso: ISSorted()
1735: @*/
1736: PetscErrorCode ISToGeneral(IS is)
1737: {
1742: if (is->ops->togeneral) {
1743: (*is->ops->togeneral)(is);
1744: } else SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Not written for this type %s",((PetscObject)is)->type_name);
1745: return(0);
1746: }
1748: /*@
1749: ISSorted - Checks the indices to determine whether they have been sorted.
1751: Collective on IS
1753: Input Parameter:
1754: . is - the index set
1756: Output Parameter:
1757: . flg - output flag, either PETSC_TRUE if the index set is sorted,
1758: or PETSC_FALSE otherwise.
1760: Notes:
1761: For parallel IS objects this only indicates if the local part of the IS
1762: is sorted. So some processors may return PETSC_TRUE while others may
1763: return PETSC_FALSE.
1765: Level: intermediate
1767: .seealso: ISSort(), ISSortRemoveDups()
1768: @*/
1769: PetscErrorCode ISSorted(IS is,PetscBool *flg)
1770: {
1776: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,flg);
1777: return(0);
1778: }
1780: /*@
1781: ISDuplicate - Creates a duplicate copy of an index set.
1783: Collective on IS
1785: Input Parmeters:
1786: . is - the index set
1788: Output Parameters:
1789: . isnew - the copy of the index set
1791: Level: beginner
1793: .seealso: ISCreateGeneral(), ISCopy()
1794: @*/
1795: PetscErrorCode ISDuplicate(IS is,IS *newIS)
1796: {
1802: (*is->ops->duplicate)(is,newIS);
1803: ISCopyInfo(is,*newIS);
1804: return(0);
1805: }
1807: /*@
1808: ISCopy - Copies an index set.
1810: Collective on IS
1812: Input Parmeters:
1813: . is - the index set
1815: Output Parameters:
1816: . isy - the copy of the index set
1818: Level: beginner
1820: .seealso: ISDuplicate()
1821: @*/
1822: PetscErrorCode ISCopy(IS is,IS isy)
1823: {
1830: if (is == isy) return(0);
1831: ISCopyInfo(is,isy);
1832: isy->max = is->max;
1833: isy->min = is->min;
1834: (*is->ops->copy)(is,isy);
1835: return(0);
1836: }
1838: /*@
1839: ISOnComm - Split a parallel IS on subcomms (usually self) or concatenate index sets on subcomms into a parallel index set
1841: Collective on IS
1843: Input Arguments:
1844: + is - index set
1845: . comm - communicator for new index set
1846: - mode - copy semantics, PETSC_USE_POINTER for no-copy if possible, otherwise PETSC_COPY_VALUES
1848: Output Arguments:
1849: . newis - new IS on comm
1851: Level: advanced
1853: Notes:
1854: It is usually desirable to create a parallel IS and look at the local part when necessary.
1856: This function is useful if serial ISs must be created independently, or to view many
1857: logically independent serial ISs.
1859: The input IS must have the same type on every process.
1860: @*/
1861: PetscErrorCode ISOnComm(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
1862: {
1864: PetscMPIInt match;
1869: MPI_Comm_compare(PetscObjectComm((PetscObject)is),comm,&match);
1870: if (mode != PETSC_COPY_VALUES && (match == MPI_IDENT || match == MPI_CONGRUENT)) {
1871: PetscObjectReference((PetscObject)is);
1872: *newis = is;
1873: } else {
1874: (*is->ops->oncomm)(is,comm,mode,newis);
1875: }
1876: return(0);
1877: }
1879: /*@
1880: ISSetBlockSize - informs an index set that it has a given block size
1882: Logicall Collective on IS
1884: Input Arguments:
1885: + is - index set
1886: - bs - block size
1888: Level: intermediate
1890: Notes:
1891: This is much like the block size for Vecs. It indicates that one can think of the indices as
1892: being in a collection of equal size blocks. For ISBlock() these collections of blocks are all contiquous
1893: within a block but this is not the case for other IS.
1894: ISBlockGetIndices() only works for ISBlock IS, not others.
1896: .seealso: ISGetBlockSize(), ISCreateBlock(), ISBlockGetIndices(),
1897: @*/
1898: PetscErrorCode ISSetBlockSize(IS is,PetscInt bs)
1899: {
1905: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_OUTOFRANGE,"Block size %D, must be positive",bs);
1906: (*is->ops->setblocksize)(is,bs);
1907: return(0);
1908: }
1910: /*@
1911: ISGetBlockSize - Returns the number of elements in a block.
1913: Not Collective
1915: Input Parameter:
1916: . is - the index set
1918: Output Parameter:
1919: . size - the number of elements in a block
1921: Level: intermediate
1923: Notes:
1924: This is much like the block size for Vecs. It indicates that one can think of the indices as
1925: being in a collection of equal size blocks. For ISBlock() these collections of blocks are all contiquous
1926: within a block but this is not the case for other IS.
1927: ISBlockGetIndices() only works for ISBlock IS, not others.
1929: .seealso: ISBlockGetSize(), ISGetSize(), ISCreateBlock(), ISSetBlockSize()
1930: @*/
1931: PetscErrorCode ISGetBlockSize(IS is,PetscInt *size)
1932: {
1936: PetscLayoutGetBlockSize(is->map, size);
1937: return(0);
1938: }
1940: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[])
1941: {
1943: PetscInt len,i;
1944: const PetscInt *ptr;
1947: ISGetLocalSize(is,&len);
1948: ISGetIndices(is,&ptr);
1949: for (i=0; i<len; i++) idx[i] = ptr[i];
1950: ISRestoreIndices(is,&ptr);
1951: return(0);
1952: }
1954: /*MC
1955: ISGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1956: The users should call ISRestoreIndicesF90() after having looked at the
1957: indices. The user should NOT change the indices.
1959: Synopsis:
1960: ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1962: Not collective
1964: Input Parameter:
1965: . x - index set
1967: Output Parameters:
1968: + xx_v - the Fortran90 pointer to the array
1969: - ierr - error code
1971: Example of Usage:
1972: .vb
1973: PetscInt, pointer xx_v(:)
1974: ....
1975: call ISGetIndicesF90(x,xx_v,ierr)
1976: a = xx_v(3)
1977: call ISRestoreIndicesF90(x,xx_v,ierr)
1978: .ve
1980: Level: intermediate
1982: .seealso: ISRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices()
1985: M*/
1987: /*MC
1988: ISRestoreIndicesF90 - Restores an index set to a usable state after
1989: a call to ISGetIndicesF90().
1991: Synopsis:
1992: ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1994: Not collective
1996: Input Parameters:
1997: + x - index set
1998: - xx_v - the Fortran90 pointer to the array
2000: Output Parameter:
2001: . ierr - error code
2004: Example of Usage:
2005: .vb
2006: PetscInt, pointer xx_v(:)
2007: ....
2008: call ISGetIndicesF90(x,xx_v,ierr)
2009: a = xx_v(3)
2010: call ISRestoreIndicesF90(x,xx_v,ierr)
2011: .ve
2013: Level: intermediate
2015: .seealso: ISGetIndicesF90(), ISGetIndices(), ISRestoreIndices()
2017: M*/
2019: /*MC
2020: ISBlockGetIndicesF90 - Accesses the elements of an index set from Fortran90.
2021: The users should call ISBlockRestoreIndicesF90() after having looked at the
2022: indices. The user should NOT change the indices.
2024: Synopsis:
2025: ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
2027: Not collective
2029: Input Parameter:
2030: . x - index set
2032: Output Parameters:
2033: + xx_v - the Fortran90 pointer to the array
2034: - ierr - error code
2035: Example of Usage:
2036: .vb
2037: PetscInt, pointer xx_v(:)
2038: ....
2039: call ISBlockGetIndicesF90(x,xx_v,ierr)
2040: a = xx_v(3)
2041: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
2042: .ve
2044: Level: intermediate
2046: .seealso: ISBlockRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices(),
2047: ISRestoreIndices()
2049: M*/
2051: /*MC
2052: ISBlockRestoreIndicesF90 - Restores an index set to a usable state after
2053: a call to ISBlockGetIndicesF90().
2055: Synopsis:
2056: ISBlockRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
2058: Not Collective
2060: Input Parameters:
2061: + x - index set
2062: - xx_v - the Fortran90 pointer to the array
2064: Output Parameter:
2065: . ierr - error code
2067: Example of Usage:
2068: .vb
2069: PetscInt, pointer xx_v(:)
2070: ....
2071: call ISBlockGetIndicesF90(x,xx_v,ierr)
2072: a = xx_v(3)
2073: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
2074: .ve
2076: Notes:
2077: Not yet supported for all F90 compilers
2079: Level: intermediate
2081: .seealso: ISBlockGetIndicesF90(), ISGetIndices(), ISRestoreIndices(), ISRestoreIndicesF90()
2083: M*/