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: }
175: /*@
176: ISCreateSubIS - Create a sub index set from a global index set selecting some components.
178: Collective on IS
180: Input Parmeters:
181: + is - the index set
182: - comps - which components we will extract from is
184: Output Parameters:
185: . subis - the new sub index set
187: Level: intermediate
189: Example usage:
190: We have an index set (is) living on 3 processes with the following values:
191: | 4 9 0 | 2 6 7 | 10 11 1|
192: and another index set (comps) used to indicate which components of is we want to take,
193: | 7 5 | 1 2 | 0 4|
194: The output index set (subis) should look like:
195: | 11 7 | 9 0 | 4 6|
197: .seealso: VecGetSubVector(), MatCreateSubMatrix()
198: @*/
199: PetscErrorCode ISCreateSubIS(IS is,IS comps,IS *subis)
200: {
201: PetscSF sf;
202: const PetscInt *is_indices,*comps_indices;
203: PetscInt *subis_indices,nroots,nleaves,*mine,i,lidx;
204: PetscMPIInt owner;
205: PetscSFNode *remote;
206: PetscErrorCode ierr;
207: MPI_Comm comm;
214: PetscObjectGetComm((PetscObject)is, &comm);
215: ISGetLocalSize(comps,&nleaves);
216: ISGetLocalSize(is,&nroots);
217: PetscMalloc1(nleaves,&remote);
218: PetscMalloc1(nleaves,&mine);
219: ISGetIndices(comps,&comps_indices);
220: /*
221: * Construct a PetscSF in which "is" data serves as roots and "subis" is leaves.
222: * Root data are sent to leaves using PetscSFBcast().
223: * */
224: for (i=0; i<nleaves; i++) {
225: mine[i] = i;
226: /* Connect a remote root with the current leaf. The value on the remote root
227: * will be received by the current local leaf.
228: * */
229: owner = -1;
230: lidx = -1;
231: PetscLayoutFindOwnerIndex(is->map,comps_indices[i],&owner,&lidx);
232: remote[i].rank = owner;
233: remote[i].index = lidx;
234: }
235: ISRestoreIndices(comps,&comps_indices);
236: PetscSFCreate(comm,&sf);
237: PetscSFSetFromOptions(sf);\
238: PetscSFSetGraph(sf,nroots,nleaves,mine,PETSC_OWN_POINTER,remote,PETSC_OWN_POINTER);
240: PetscMalloc1(nleaves,&subis_indices);
241: ISGetIndices(is, &is_indices);
242: PetscSFBcastBegin(sf,MPIU_INT,is_indices,subis_indices,MPI_REPLACE);
243: PetscSFBcastEnd(sf,MPIU_INT,is_indices,subis_indices,MPI_REPLACE);
244: ISRestoreIndices(is,&is_indices);
245: PetscSFDestroy(&sf);
246: ISCreateGeneral(comm,nleaves,subis_indices,PETSC_OWN_POINTER,subis);
247: return(0);
248: }
250: /*@
251: ISClearInfoCache - clear the cache of computed index set properties
253: Not collective
255: Input Parameters:
256: + is - the index set
257: - clear_permanent_local - whether to remove the permanent status of local properties
259: NOTE: because all processes must agree on the global permanent status of a property,
260: the permanent status can only be changed with ISSetInfo(), because this routine is not collective
262: Level: developer
264: .seealso: ISInfo, ISInfoType, ISSetInfo(), ISClearInfoCache()
266: @*/
267: PetscErrorCode ISClearInfoCache(IS is, PetscBool clear_permanent_local)
268: {
269: PetscInt i, j;
274: for (i = 0; i < IS_INFO_MAX; i++) {
275: if (clear_permanent_local) is->info_permanent[IS_LOCAL][i] = PETSC_FALSE;
276: for (j = 0; j < 2; j++) {
277: if (!is->info_permanent[j][i]) is->info[j][i] = IS_INFO_UNKNOWN;
278: }
279: }
280: return(0);
281: }
283: static PetscErrorCode ISSetInfo_Internal(IS is, ISInfo info, ISInfoType type, ISInfoBool ipermanent, PetscBool flg)
284: {
285: ISInfoBool iflg = flg ? IS_INFO_TRUE : IS_INFO_FALSE;
286: PetscInt itype = (type == IS_LOCAL) ? 0 : 1;
287: PetscBool permanent_set = (ipermanent == IS_INFO_UNKNOWN) ? PETSC_FALSE : PETSC_TRUE;
288: PetscBool permanent = (ipermanent == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
291: /* set this property */
292: is->info[itype][(int)info] = iflg;
293: if (permanent_set) is->info_permanent[itype][(int)info] = permanent;
294: /* set implications */
295: switch (info) {
296: case IS_SORTED:
297: if (flg && type == IS_GLOBAL) { /* an array that is globally sorted is also locally sorted */
298: is->info[IS_LOCAL][(int)info] = IS_INFO_TRUE;
299: /* global permanence implies local permanence */
300: if (permanent_set && permanent) is->info_permanent[IS_LOCAL][(int)info] = PETSC_TRUE;
301: }
302: if (!flg) { /* if an array is not sorted, it cannot be an interval or the identity */
303: is->info[itype][IS_INTERVAL] = IS_INFO_FALSE;
304: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
305: if (permanent_set) {
306: is->info_permanent[itype][IS_INTERVAL] = permanent;
307: is->info_permanent[itype][IS_IDENTITY] = permanent;
308: }
309: }
310: break;
311: case IS_UNIQUE:
312: if (flg && type == IS_GLOBAL) { /* an array that is globally unique is also locally unique */
313: is->info[IS_LOCAL][(int)info] = IS_INFO_TRUE;
314: /* global permanence implies local permanence */
315: if (permanent_set && permanent) is->info_permanent[IS_LOCAL][(int)info] = PETSC_TRUE;
316: }
317: if (!flg) { /* if an array is not unique, it cannot be a permutation, and interval, or the identity */
318: is->info[itype][IS_PERMUTATION] = IS_INFO_FALSE;
319: is->info[itype][IS_INTERVAL] = IS_INFO_FALSE;
320: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
321: if (permanent_set) {
322: is->info_permanent[itype][IS_PERMUTATION] = permanent;
323: is->info_permanent[itype][IS_INTERVAL] = permanent;
324: is->info_permanent[itype][IS_IDENTITY] = permanent;
325: }
326: }
327: break;
328: case IS_PERMUTATION:
329: if (flg) { /* an array that is a permutation is unique and is unique locally */
330: is->info[itype][IS_UNIQUE] = IS_INFO_TRUE;
331: is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
332: if (permanent_set && permanent) {
333: is->info_permanent[itype][IS_UNIQUE] = PETSC_TRUE;
334: is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
335: }
336: } else { /* an array that is not a permutation cannot be the identity */
337: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
338: if (permanent_set) is->info_permanent[itype][IS_IDENTITY] = permanent;
339: }
340: break;
341: case IS_INTERVAL:
342: if (flg) { /* an array that is an interval is sorted and unique */
343: is->info[itype][IS_SORTED] = IS_INFO_TRUE;
344: is->info[IS_LOCAL][IS_SORTED] = IS_INFO_TRUE;
345: is->info[itype][IS_UNIQUE] = IS_INFO_TRUE;
346: is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
347: if (permanent_set && permanent) {
348: is->info_permanent[itype][IS_SORTED] = PETSC_TRUE;
349: is->info_permanent[IS_LOCAL][IS_SORTED] = PETSC_TRUE;
350: is->info_permanent[itype][IS_UNIQUE] = PETSC_TRUE;
351: is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
352: }
353: } else { /* an array that is not an interval cannot be the identity */
354: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
355: if (permanent_set) is->info_permanent[itype][IS_IDENTITY] = permanent;
356: }
357: break;
358: case IS_IDENTITY:
359: if (flg) { /* an array that is the identity is sorted, unique, an interval, and a permutation */
360: is->info[itype][IS_SORTED] = IS_INFO_TRUE;
361: is->info[IS_LOCAL][IS_SORTED] = IS_INFO_TRUE;
362: is->info[itype][IS_UNIQUE] = IS_INFO_TRUE;
363: is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
364: is->info[itype][IS_PERMUTATION] = IS_INFO_TRUE;
365: is->info[itype][IS_INTERVAL] = IS_INFO_TRUE;
366: is->info[IS_LOCAL][IS_INTERVAL] = IS_INFO_TRUE;
367: if (permanent_set && permanent) {
368: is->info_permanent[itype][IS_SORTED] = PETSC_TRUE;
369: is->info_permanent[IS_LOCAL][IS_SORTED] = PETSC_TRUE;
370: is->info_permanent[itype][IS_UNIQUE] = PETSC_TRUE;
371: is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
372: is->info_permanent[itype][IS_PERMUTATION] = PETSC_TRUE;
373: is->info_permanent[itype][IS_INTERVAL] = PETSC_TRUE;
374: is->info_permanent[IS_LOCAL][IS_INTERVAL] = PETSC_TRUE;
375: }
376: }
377: break;
378: default:
379: if (type == IS_LOCAL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
380: else SETERRQ(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
381: }
382: return(0);
383: }
385: /*@
386: ISSetInfo - Set known information about an index set.
388: Logically Collective on IS if type is IS_GLOBAL
390: Input Parameters:
391: + is - the index set
392: . info - describing a property of the index set, one of those listed below,
393: . type - IS_LOCAL if the information describes the local portion of the index set,
394: IS_GLOBAL if it describes the whole index set
395: . permanent - PETSC_TRUE if it is known that the property will persist through changes to the index set, PETSC_FALSE otherwise
396: If the user sets a property as permanently known, it will bypass computation of that property
397: - flg - set the described property as true (PETSC_TRUE) or false (PETSC_FALSE)
399: Info Describing IS Structure:
400: + IS_SORTED - the [local part of the] index set is sorted in ascending order
401: . IS_UNIQUE - each entry in the [local part of the] index set is unique
402: . 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
403: . IS_INTERVAL - the [local part of the] index set is equal to a contiguous range of integers {f, f + 1, ..., f + N-1}
404: - IS_IDENTITY - the [local part of the] index set is equal to the integers {0, 1, ..., N-1}
406: Notes:
407: If type is IS_GLOBAL, all processes that share the index set must pass the same value in flg
409: It is possible to set a property with ISSetInfo() that contradicts what would be previously computed with ISGetInfo()
411: Level: advanced
413: .seealso: ISInfo, ISInfoType, IS
415: @*/
416: PetscErrorCode ISSetInfo(IS is, ISInfo info, ISInfoType type, PetscBool permanent, PetscBool flg)
417: {
418: MPI_Comm comm, errcomm;
419: PetscMPIInt size;
425: comm = PetscObjectComm((PetscObject)is);
426: if (type == IS_GLOBAL) {
430: errcomm = comm;
431: } else {
432: errcomm = PETSC_COMM_SELF;
433: }
435: 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);
437: MPI_Comm_size(comm, &size);
438: /* do not use global values if size == 1: it makes it easier to keep the implications straight */
439: if (size == 1) type = IS_LOCAL;
440: ISSetInfo_Internal(is, info, type, permanent ? IS_INFO_TRUE : IS_INFO_FALSE, flg);
441: return(0);
442: }
444: static PetscErrorCode ISGetInfo_Sorted(IS is, ISInfoType type, PetscBool *flg)
445: {
446: MPI_Comm comm;
447: PetscMPIInt size, rank;
451: comm = PetscObjectComm((PetscObject)is);
452: MPI_Comm_size(comm, &size);
453: MPI_Comm_size(comm, &rank);
454: if (type == IS_GLOBAL && is->ops->sortedglobal) {
455: (*is->ops->sortedglobal)(is,flg);
456: } else {
457: PetscBool sortedLocal = PETSC_FALSE;
459: /* determine if the array is locally sorted */
460: if (type == IS_GLOBAL && size > 1) {
461: /* call ISGetInfo so that a cached value will be used if possible */
462: ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sortedLocal);
463: } else if (is->ops->sortedlocal) {
464: (*is->ops->sortedlocal)(is,&sortedLocal);
465: } else {
466: /* default: get the local indices and directly check */
467: const PetscInt *idx;
468: PetscInt n;
470: ISGetIndices(is, &idx);
471: ISGetLocalSize(is, &n);
472: PetscSortedInt(n, idx, &sortedLocal);
473: ISRestoreIndices(is, &idx);
474: }
476: if (type == IS_LOCAL || size == 1) {
477: *flg = sortedLocal;
478: } else {
479: MPI_Allreduce(&sortedLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
480: if (*flg) {
481: PetscInt n, min = PETSC_MAX_INT, max = PETSC_MIN_INT;
482: PetscInt maxprev;
484: ISGetLocalSize(is, &n);
485: if (n) {ISGetMinMax(is, &min, &max);}
486: maxprev = PETSC_MIN_INT;
487: MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm);
488: if (rank && (maxprev > min)) sortedLocal = PETSC_FALSE;
489: MPI_Allreduce(&sortedLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
490: }
491: }
492: }
493: return(0);
494: }
496: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[]);
498: static PetscErrorCode ISGetInfo_Unique(IS is, ISInfoType type, PetscBool *flg)
499: {
500: MPI_Comm comm;
501: PetscMPIInt size, rank;
502: PetscInt i;
506: comm = PetscObjectComm((PetscObject)is);
507: MPI_Comm_size(comm, &size);
508: MPI_Comm_size(comm, &rank);
509: if (type == IS_GLOBAL && is->ops->uniqueglobal) {
510: (*is->ops->uniqueglobal)(is,flg);
511: } else {
512: PetscBool uniqueLocal;
513: PetscInt n = -1;
514: PetscInt *idx = NULL;
516: /* determine if the array is locally unique */
517: if (type == IS_GLOBAL && size > 1) {
518: /* call ISGetInfo so that a cached value will be used if possible */
519: ISGetInfo(is, IS_UNIQUE, IS_LOCAL, PETSC_TRUE, &uniqueLocal);
520: } else if (is->ops->uniquelocal) {
521: (*is->ops->uniquelocal)(is,&uniqueLocal);
522: } else {
523: /* default: get the local indices and directly check */
524: uniqueLocal = PETSC_TRUE;
525: ISGetLocalSize(is, &n);
526: PetscMalloc1(n, &idx);
527: ISGetIndicesCopy(is, idx);
528: PetscIntSortSemiOrdered(n, idx);
529: for (i = 1; i < n; i++) if (idx[i] == idx[i-1]) break;
530: if (i < n) uniqueLocal = PETSC_FALSE;
531: }
533: PetscFree(idx);
534: if (type == IS_LOCAL || size == 1) {
535: *flg = uniqueLocal;
536: } else {
537: MPI_Allreduce(&uniqueLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
538: if (*flg) {
539: PetscInt min = PETSC_MAX_INT, max = PETSC_MIN_INT, maxprev;
541: if (!idx) {
542: ISGetLocalSize(is, &n);
543: PetscMalloc1(n, &idx);
544: ISGetIndicesCopy(is, idx);
545: }
546: PetscParallelSortInt(is->map, is->map, idx, idx);
547: if (n) {
548: min = idx[0];
549: max = idx[n - 1];
550: }
551: for (i = 1; i < n; i++) if (idx[i] == idx[i-1]) break;
552: if (i < n) uniqueLocal = PETSC_FALSE;
553: maxprev = PETSC_MIN_INT;
554: MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm);
555: if (rank && (maxprev == min)) uniqueLocal = PETSC_FALSE;
556: MPI_Allreduce(&uniqueLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
557: }
558: }
559: PetscFree(idx);
560: }
561: return(0);
562: }
564: static PetscErrorCode ISGetInfo_Permutation(IS is, ISInfoType type, PetscBool *flg)
565: {
566: MPI_Comm comm;
567: PetscMPIInt size, rank;
571: comm = PetscObjectComm((PetscObject)is);
572: MPI_Comm_size(comm, &size);
573: MPI_Comm_size(comm, &rank);
574: if (type == IS_GLOBAL && is->ops->permglobal) {
575: (*is->ops->permglobal)(is,flg);
576: } else if (type == IS_LOCAL && is->ops->permlocal) {
577: (*is->ops->permlocal)(is,flg);
578: } else {
579: PetscBool permLocal;
580: PetscInt n, i, rStart;
581: PetscInt *idx;
583: ISGetLocalSize(is, &n);
584: PetscMalloc1(n, &idx);
585: ISGetIndicesCopy(is, idx);
586: if (type == IS_GLOBAL) {
587: PetscParallelSortInt(is->map, is->map, idx, idx);
588: PetscLayoutGetRange(is->map, &rStart, NULL);
589: } else {
590: PetscIntSortSemiOrdered(n, idx);
591: rStart = 0;
592: }
593: permLocal = PETSC_TRUE;
594: for (i = 0; i < n; i++) {
595: if (idx[i] != rStart + i) break;
596: }
597: if (i < n) permLocal = PETSC_FALSE;
598: if (type == IS_LOCAL || size == 1) {
599: *flg = permLocal;
600: } else {
601: MPI_Allreduce(&permLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
602: }
603: PetscFree(idx);
604: }
605: return(0);
606: }
608: static PetscErrorCode ISGetInfo_Interval(IS is, ISInfoType type, PetscBool *flg)
609: {
610: MPI_Comm comm;
611: PetscMPIInt size, rank;
612: PetscInt i;
616: comm = PetscObjectComm((PetscObject)is);
617: MPI_Comm_size(comm, &size);
618: MPI_Comm_size(comm, &rank);
619: if (type == IS_GLOBAL && is->ops->intervalglobal) {
620: (*is->ops->intervalglobal)(is,flg);
621: } else {
622: PetscBool intervalLocal;
624: /* determine if the array is locally an interval */
625: if (type == IS_GLOBAL && size > 1) {
626: /* call ISGetInfo so that a cached value will be used if possible */
627: ISGetInfo(is, IS_INTERVAL, IS_LOCAL, PETSC_TRUE, &intervalLocal);
628: } else if (is->ops->intervallocal) {
629: (*is->ops->intervallocal)(is,&intervalLocal);
630: } else {
631: PetscInt n;
632: const PetscInt *idx;
633: /* default: get the local indices and directly check */
634: intervalLocal = PETSC_TRUE;
635: ISGetLocalSize(is, &n);
636: PetscMalloc1(n, &idx);
637: ISGetIndices(is, &idx);
638: for (i = 1; i < n; i++) if (idx[i] != idx[i-1] + 1) break;
639: if (i < n) intervalLocal = PETSC_FALSE;
640: ISRestoreIndices(is, &idx);
641: }
643: if (type == IS_LOCAL || size == 1) {
644: *flg = intervalLocal;
645: } else {
646: MPI_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
647: if (*flg) {
648: PetscInt n, min = PETSC_MAX_INT, max = PETSC_MIN_INT;
649: PetscInt maxprev;
651: ISGetLocalSize(is, &n);
652: if (n) {ISGetMinMax(is, &min, &max);}
653: maxprev = PETSC_MIN_INT;
654: MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm);
655: if (rank && n && (maxprev != min - 1)) intervalLocal = PETSC_FALSE;
656: MPI_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
657: }
658: }
659: }
660: return(0);
661: }
663: static PetscErrorCode ISGetInfo_Identity(IS is, ISInfoType type, PetscBool *flg)
664: {
665: MPI_Comm comm;
666: PetscMPIInt size, rank;
670: comm = PetscObjectComm((PetscObject)is);
671: MPI_Comm_size(comm, &size);
672: MPI_Comm_size(comm, &rank);
673: if (type == IS_GLOBAL && is->ops->intervalglobal) {
674: PetscBool isinterval;
676: (*is->ops->intervalglobal)(is,&isinterval);
677: *flg = PETSC_FALSE;
678: if (isinterval) {
679: PetscInt min;
681: ISGetMinMax(is, &min, NULL);
682: MPI_Bcast(&min, 1, MPIU_INT, 0, comm);
683: if (min == 0) *flg = PETSC_TRUE;
684: }
685: } else if (type == IS_LOCAL && is->ops->intervallocal) {
686: PetscBool isinterval;
688: (*is->ops->intervallocal)(is,&isinterval);
689: *flg = PETSC_FALSE;
690: if (isinterval) {
691: PetscInt min;
693: ISGetMinMax(is, &min, NULL);
694: if (min == 0) *flg = PETSC_TRUE;
695: }
696: } else {
697: PetscBool identLocal;
698: PetscInt n, i, rStart;
699: const PetscInt *idx;
701: ISGetLocalSize(is, &n);
702: ISGetIndices(is, &idx);
703: PetscLayoutGetRange(is->map, &rStart, NULL);
704: identLocal = PETSC_TRUE;
705: for (i = 0; i < n; i++) {
706: if (idx[i] != rStart + i) break;
707: }
708: if (i < n) identLocal = PETSC_FALSE;
709: if (type == IS_LOCAL || size == 1) {
710: *flg = identLocal;
711: } else {
712: MPI_Allreduce(&identLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
713: }
714: ISRestoreIndices(is, &idx);
715: }
716: return(0);
717: }
719: /*@
720: ISGetInfo - Determine whether an index set satisfies a given property
722: 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())
724: Input Parameters:
725: + is - the index set
726: . info - describing a property of the index set, one of those listed in the documentation of ISSetInfo()
727: . 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
728: - type - whether the property is local (IS_LOCAL) or global (IS_GLOBAL)
730: Output Parameter:
731: . flg - wheter the property is true (PETSC_TRUE) or false (PETSC_FALSE)
733: 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().
735: Level: advanced
737: .seealso: ISInfo, ISInfoType, ISSetInfo(), ISClearInfoCache()
739: @*/
740: PetscErrorCode ISGetInfo(IS is, ISInfo info, ISInfoType type, PetscBool compute, PetscBool *flg)
741: {
742: MPI_Comm comm, errcomm;
743: PetscMPIInt rank, size;
744: PetscInt itype;
745: PetscBool hasprop;
746: PetscBool infer;
752: comm = PetscObjectComm((PetscObject)is);
753: if (type == IS_GLOBAL) {
755: errcomm = comm;
756: } else {
757: errcomm = PETSC_COMM_SELF;
758: }
760: MPI_Comm_size(comm, &size);
761: MPI_Comm_rank(comm, &rank);
763: 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);
764: if (size == 1) type = IS_LOCAL;
765: itype = (type == IS_LOCAL) ? 0 : 1;
766: hasprop = PETSC_FALSE;
767: infer = PETSC_FALSE;
768: if (is->info_permanent[itype][(int)info]) {
769: hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
770: infer = PETSC_TRUE;
771: } else if ((itype == IS_LOCAL) && (is->info[IS_LOCAL][info] != IS_INFO_UNKNOWN)) {
772: /* we can cache local properties as long as we clear them when the IS changes */
773: /* NOTE: we only cache local values because there is no ISAssemblyBegin()/ISAssemblyEnd(),
774: so we have no way of knowing when a cached value has been invalidated by changes on a different process */
775: hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
776: infer = PETSC_TRUE;
777: } else if (compute) {
778: switch (info) {
779: case IS_SORTED:
780: ISGetInfo_Sorted(is, type, &hasprop);
781: break;
782: case IS_UNIQUE:
783: ISGetInfo_Unique(is, type, &hasprop);
784: break;
785: case IS_PERMUTATION:
786: ISGetInfo_Permutation(is, type, &hasprop);
787: break;
788: case IS_INTERVAL:
789: ISGetInfo_Interval(is, type, &hasprop);
790: break;
791: case IS_IDENTITY:
792: ISGetInfo_Identity(is, type, &hasprop);
793: break;
794: default:
795: SETERRQ(errcomm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
796: }
797: infer = PETSC_TRUE;
798: }
799: /* call ISSetInfo_Internal to keep all of the implications straight */
800: if (infer) {ISSetInfo_Internal(is, info, type, IS_INFO_UNKNOWN, hasprop);}
801: *flg = hasprop;
802: return(0);
803: }
805: static PetscErrorCode ISCopyInfo(IS source, IS dest)
806: {
810: PetscArraycpy(&dest->info[0], &source->info[0], 2);
811: PetscArraycpy(&dest->info_permanent[0], &source->info_permanent[0], 2);
812: return(0);
813: }
815: /*@
816: ISIdentity - Determines whether index set is the identity mapping.
818: Collective on IS
820: Input Parmeters:
821: . is - the index set
823: Output Parameters:
824: . ident - PETSC_TRUE if an identity, else PETSC_FALSE
826: Level: intermediate
828: Note: If ISSetIdentity() (or ISSetInfo() for a permanent property) has been called,
829: ISIdentity() will return its answer without communication between processes, but
830: otherwise the output ident will be computed from ISGetInfo(),
831: which may require synchronization on the communicator of IS. To avoid this computation,
832: call ISGetInfo() directly with the compute flag set to PETSC_FALSE, and ident will be assumed false.
834: .seealso: ISSetIdentity(), ISGetInfo()
835: @*/
836: PetscErrorCode ISIdentity(IS is,PetscBool *ident)
837: {
843: ISGetInfo(is,IS_IDENTITY,IS_GLOBAL,PETSC_TRUE,ident);
844: return(0);
845: }
847: /*@
848: ISSetIdentity - Informs the index set that it is an identity.
850: Logically Collective on IS
852: Input Parmeters:
853: . is - the index set
855: Level: intermediate
857: Note: The IS will be considered the identity permanently, even if indices have been changes (for example, with
858: ISGeneralSetIndices()). It's a good idea to only set this property if the IS will not change in the future.
859: To clear this property, use ISClearInfoCache().
861: .seealso: ISIdentity(), ISSetInfo(), ISClearInfoCache()
862: @*/
863: PetscErrorCode ISSetIdentity(IS is)
864: {
869: ISSetInfo(is,IS_IDENTITY,IS_GLOBAL,PETSC_TRUE,PETSC_TRUE);
870: return(0);
871: }
873: /*@
874: ISContiguousLocal - Locates an index set with contiguous range within a global range, if possible
876: Not Collective
878: Input Parmeters:
879: + is - the index set
880: . gstart - global start
881: - gend - global end
883: Output Parameters:
884: + start - start of contiguous block, as an offset from gstart
885: - contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
887: Level: developer
889: .seealso: ISGetLocalSize(), VecGetOwnershipRange()
890: @*/
891: PetscErrorCode ISContiguousLocal(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
892: {
899: *start = -1;
900: *contig = PETSC_FALSE;
901: if (is->ops->contiguous) {
902: (*is->ops->contiguous)(is,gstart,gend,start,contig);
903: }
904: return(0);
905: }
907: /*@
908: ISPermutation - PETSC_TRUE or PETSC_FALSE depending on whether the
909: index set has been declared to be a permutation.
911: Logically Collective on IS
913: Input Parmeters:
914: . is - the index set
916: Output Parameters:
917: . perm - PETSC_TRUE if a permutation, else PETSC_FALSE
919: Level: intermediate
921: Note: If it is not alread known that the IS is a permutation (if ISSetPermutation()
922: or ISSetInfo() has not been called), this routine will not attempt to compute
923: whether the index set is a permutation and will assume perm is PETSC_FALSE.
924: To compute the value when it is not already known, use ISGetInfo() with
925: the compute flag set to PETSC_TRUE.
927: .seealso: ISSetPermutation(), ISGetInfo()
928: @*/
929: PetscErrorCode ISPermutation(IS is,PetscBool *perm)
930: {
936: ISGetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_FALSE,perm);
937: return(0);
938: }
940: /*@
941: ISSetPermutation - Informs the index set that it is a permutation.
943: Logically Collective on IS
945: Input Parmeters:
946: . is - the index set
948: Level: intermediate
950: The debug version of the libraries (./configure --with-debugging=1) checks if the
951: index set is actually a permutation. The optimized version just believes you.
953: Note: The IS will be considered a permutation permanently, even if indices have been changes (for example, with
954: ISGeneralSetIndices()). It's a good idea to only set this property if the IS will not change in the future.
955: To clear this property, use ISClearInfoCache().
957: .seealso: ISPermutation(), ISSetInfo(), ISClearInfoCache().
958: @*/
959: PetscErrorCode ISSetPermutation(IS is)
960: {
965: if (PetscDefined(USE_DEBUG)) {
966: PetscMPIInt size;
968: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
969: if (size == 1) {
970: PetscInt i,n,*idx;
971: const PetscInt *iidx;
973: ISGetSize(is,&n);
974: PetscMalloc1(n,&idx);
975: ISGetIndices(is,&iidx);
976: PetscArraycpy(idx,iidx,n);
977: PetscIntSortSemiOrdered(n,idx);
978: for (i=0; i<n; i++) {
979: if (idx[i] != i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Index set is not a permutation");
980: }
981: PetscFree(idx);
982: ISRestoreIndices(is,&iidx);
983: }
984: }
985: ISSetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_TRUE,PETSC_TRUE);
986: return(0);
987: }
989: /*@C
990: ISDestroy - Destroys an index set.
992: Collective on IS
994: Input Parameters:
995: . is - the index set
997: Level: beginner
999: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlocked()
1000: @*/
1001: PetscErrorCode ISDestroy(IS *is)
1002: {
1006: if (!*is) return(0);
1008: if (--((PetscObject)(*is))->refct > 0) {*is = NULL; return(0);}
1009: if ((*is)->complement) {
1010: PetscInt refcnt;
1011: PetscObjectGetReference((PetscObject)((*is)->complement), &refcnt);
1012: if (refcnt > 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Nonlocal IS has not been restored");
1013: ISDestroy(&(*is)->complement);
1014: }
1015: if ((*is)->ops->destroy) {
1016: (*(*is)->ops->destroy)(*is);
1017: }
1018: PetscLayoutDestroy(&(*is)->map);
1019: /* Destroy local representations of offproc data. */
1020: PetscFree((*is)->total);
1021: PetscFree((*is)->nonlocal);
1022: PetscHeaderDestroy(is);
1023: return(0);
1024: }
1026: /*@
1027: ISInvertPermutation - Creates a new permutation that is the inverse of
1028: a given permutation.
1030: Collective on IS
1032: Input Parameter:
1033: + is - the index set
1034: - nlocal - number of indices on this processor in result (ignored for 1 proccessor) or
1035: use PETSC_DECIDE
1037: Output Parameter:
1038: . isout - the inverse permutation
1040: Level: intermediate
1042: Notes:
1043: For parallel index sets this does the complete parallel permutation, but the
1044: code is not efficient for huge index sets (10,000,000 indices).
1046: @*/
1047: PetscErrorCode ISInvertPermutation(IS is,PetscInt nlocal,IS *isout)
1048: {
1049: PetscBool isperm, isidentity, issame;
1055: ISGetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_TRUE,&isperm);
1056: if (!isperm) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_WRONG,"Not a permutation");
1057: ISGetInfo(is,IS_IDENTITY,IS_GLOBAL,PETSC_TRUE,&isidentity);
1058: issame = PETSC_FALSE;
1059: if (isidentity) {
1060: PetscInt n;
1061: PetscBool isallsame;
1063: ISGetLocalSize(is, &n);
1064: issame = (PetscBool) (n == nlocal);
1065: MPI_Allreduce(&issame, &isallsame, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is));
1066: issame = isallsame;
1067: }
1068: if (issame) {
1069: ISDuplicate(is,isout);
1070: } else {
1071: (*is->ops->invertpermutation)(is,nlocal,isout);
1072: ISSetPermutation(*isout);
1073: }
1074: return(0);
1075: }
1077: /*@
1078: ISGetSize - Returns the global length of an index set.
1080: Not Collective
1082: Input Parameter:
1083: . is - the index set
1085: Output Parameter:
1086: . size - the global size
1088: Level: beginner
1090: @*/
1091: PetscErrorCode ISGetSize(IS is,PetscInt *size)
1092: {
1096: *size = is->map->N;
1097: return(0);
1098: }
1100: /*@
1101: ISGetLocalSize - Returns the local (processor) length of an index set.
1103: Not Collective
1105: Input Parameter:
1106: . is - the index set
1108: Output Parameter:
1109: . size - the local size
1111: Level: beginner
1113: @*/
1114: PetscErrorCode ISGetLocalSize(IS is,PetscInt *size)
1115: {
1119: *size = is->map->n;
1120: return(0);
1121: }
1123: /*@
1124: ISGetLayout - get PetscLayout describing index set layout
1126: Not Collective
1128: Input Arguments:
1129: . is - the index set
1131: Output Arguments:
1132: . map - the layout
1134: Level: developer
1136: .seealso: ISGetSize(), ISGetLocalSize()
1137: @*/
1138: PetscErrorCode ISGetLayout(IS is,PetscLayout *map)
1139: {
1144: *map = is->map;
1145: return(0);
1146: }
1148: /*@C
1149: ISGetIndices - Returns a pointer to the indices. The user should call
1150: ISRestoreIndices() after having looked at the indices. The user should
1151: NOT change the indices.
1153: Not Collective
1155: Input Parameter:
1156: . is - the index set
1158: Output Parameter:
1159: . ptr - the location to put the pointer to the indices
1161: Fortran Note:
1162: This routine has two different interfaces from Fortran; the first is not recommend, it does not require Fortran 90
1163: $ IS is
1164: $ integer is_array(1)
1165: $ PetscOffset i_is
1166: $ int ierr
1167: $ call ISGetIndices(is,is_array,i_is,ierr)
1168: $
1169: $ Access first local entry in list
1170: $ value = is_array(i_is + 1)
1171: $
1172: $ ...... other code
1173: $ call ISRestoreIndices(is,is_array,i_is,ierr)
1174: The second Fortran interface is recommended.
1175: $ use petscisdef
1176: $ PetscInt, pointer :: array(:)
1177: $ PetscErrorCode ierr
1178: $ IS i
1179: $ call ISGetIndicesF90(i,array,ierr)
1181: See the Fortran chapter of the users manual and
1182: petsc/src/is/[tutorials,tests] for details.
1184: Level: intermediate
1186: .seealso: ISRestoreIndices(), ISGetIndicesF90()
1187: @*/
1188: PetscErrorCode ISGetIndices(IS is,const PetscInt *ptr[])
1189: {
1195: (*is->ops->getindices)(is,ptr);
1196: return(0);
1197: }
1199: /*@C
1200: ISGetMinMax - Gets the minimum and maximum values in an IS
1202: Not Collective
1204: Input Parameter:
1205: . is - the index set
1207: Output Parameter:
1208: + min - the minimum value
1209: - max - the maximum value
1211: Level: intermediate
1213: Notes:
1214: Empty index sets return min=PETSC_MAX_INT and max=PETSC_MIN_INT.
1215: In parallel, it returns the min and max of the local portion of the IS
1217: .seealso: ISGetIndices(), ISRestoreIndices(), ISGetIndicesF90()
1218: @*/
1219: PetscErrorCode ISGetMinMax(IS is,PetscInt *min,PetscInt *max)
1220: {
1223: if (min) *min = is->min;
1224: if (max) *max = is->max;
1225: return(0);
1226: }
1228: /*@
1229: ISLocate - determine the location of an index within the local component of an index set
1231: Not Collective
1233: Input Parameter:
1234: + is - the index set
1235: - key - the search key
1237: Output Parameter:
1238: . location - if >= 0, a location within the index set that is equal to the key, otherwise the key is not in the index set
1240: Level: intermediate
1241: @*/
1242: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
1243: {
1247: if (is->ops->locate) {
1248: (*is->ops->locate)(is,key,location);
1249: } else {
1250: PetscInt numIdx;
1251: PetscBool sorted;
1252: const PetscInt *idx;
1254: ISGetLocalSize(is,&numIdx);
1255: ISGetIndices(is,&idx);
1256: ISSorted(is,&sorted);
1257: if (sorted) {
1258: PetscFindInt(key,numIdx,idx,location);
1259: } else {
1260: PetscInt i;
1262: *location = -1;
1263: for (i = 0; i < numIdx; i++) {
1264: if (idx[i] == key) {
1265: *location = i;
1266: break;
1267: }
1268: }
1269: }
1270: ISRestoreIndices(is,&idx);
1271: }
1272: return(0);
1273: }
1275: /*@C
1276: ISRestoreIndices - Restores an index set to a usable state after a call
1277: to ISGetIndices().
1279: Not Collective
1281: Input Parameters:
1282: + is - the index set
1283: - ptr - the pointer obtained by ISGetIndices()
1285: Fortran Note:
1286: This routine is used differently from Fortran
1287: $ IS is
1288: $ integer is_array(1)
1289: $ PetscOffset i_is
1290: $ int ierr
1291: $ call ISGetIndices(is,is_array,i_is,ierr)
1292: $
1293: $ Access first local entry in list
1294: $ value = is_array(i_is + 1)
1295: $
1296: $ ...... other code
1297: $ call ISRestoreIndices(is,is_array,i_is,ierr)
1299: See the Fortran chapter of the users manual and
1300: petsc/src/vec/is/tests for details.
1302: Level: intermediate
1304: Note:
1305: This routine zeros out ptr. This is to prevent accidental us of the array after it has been restored.
1307: .seealso: ISGetIndices(), ISRestoreIndicesF90()
1308: @*/
1309: PetscErrorCode ISRestoreIndices(IS is,const PetscInt *ptr[])
1310: {
1316: if (is->ops->restoreindices) {
1317: (*is->ops->restoreindices)(is,ptr);
1318: }
1319: return(0);
1320: }
1322: static PetscErrorCode ISGatherTotal_Private(IS is)
1323: {
1325: PetscInt i,n,N;
1326: const PetscInt *lindices;
1327: MPI_Comm comm;
1328: PetscMPIInt rank,size,*sizes = NULL,*offsets = NULL,nn;
1333: PetscObjectGetComm((PetscObject)is,&comm);
1334: MPI_Comm_size(comm,&size);
1335: MPI_Comm_rank(comm,&rank);
1336: ISGetLocalSize(is,&n);
1337: PetscMalloc2(size,&sizes,size,&offsets);
1339: PetscMPIIntCast(n,&nn);
1340: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
1341: offsets[0] = 0;
1342: for (i=1; i<size; ++i) offsets[i] = offsets[i-1] + sizes[i-1];
1343: N = offsets[size-1] + sizes[size-1];
1345: PetscMalloc1(N,&(is->total));
1346: ISGetIndices(is,&lindices);
1347: MPI_Allgatherv((void*)lindices,nn,MPIU_INT,is->total,sizes,offsets,MPIU_INT,comm);
1348: ISRestoreIndices(is,&lindices);
1349: is->local_offset = offsets[rank];
1350: PetscFree2(sizes,offsets);
1351: return(0);
1352: }
1354: /*@C
1355: ISGetTotalIndices - Retrieve an array containing all indices across the communicator.
1357: Collective on IS
1359: Input Parameter:
1360: . is - the index set
1362: Output Parameter:
1363: . indices - total indices with rank 0 indices first, and so on; total array size is
1364: the same as returned with ISGetSize().
1366: Level: intermediate
1368: Notes:
1369: this is potentially nonscalable, but depends on the size of the total index set
1370: and the size of the communicator. This may be feasible for index sets defined on
1371: subcommunicators, such that the set size does not grow with PETSC_WORLD_COMM.
1372: Note also that there is no way to tell where the local part of the indices starts
1373: (use ISGetIndices() and ISGetNonlocalIndices() to retrieve just the local and just
1374: the nonlocal part (complement), respectively).
1376: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices(), ISGetSize()
1377: @*/
1378: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
1379: {
1381: PetscMPIInt size;
1386: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
1387: if (size == 1) {
1388: (*is->ops->getindices)(is,indices);
1389: } else {
1390: if (!is->total) {
1391: ISGatherTotal_Private(is);
1392: }
1393: *indices = is->total;
1394: }
1395: return(0);
1396: }
1398: /*@C
1399: ISRestoreTotalIndices - Restore the index array obtained with ISGetTotalIndices().
1401: Not Collective.
1403: Input Parameter:
1404: + is - the index set
1405: - indices - index array; must be the array obtained with ISGetTotalIndices()
1407: Level: intermediate
1409: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices()
1410: @*/
1411: PetscErrorCode ISRestoreTotalIndices(IS is, const PetscInt *indices[])
1412: {
1414: PetscMPIInt size;
1419: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
1420: if (size == 1) {
1421: (*is->ops->restoreindices)(is,indices);
1422: } else {
1423: 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.");
1424: }
1425: return(0);
1426: }
1427: /*@C
1428: ISGetNonlocalIndices - Retrieve an array of indices from remote processors
1429: in this communicator.
1431: Collective on IS
1433: Input Parameter:
1434: . is - the index set
1436: Output Parameter:
1437: . indices - indices with rank 0 indices first, and so on, omitting
1438: the current rank. Total number of indices is the difference
1439: total and local, obtained with ISGetSize() and ISGetLocalSize(),
1440: respectively.
1442: Level: intermediate
1444: Notes:
1445: restore the indices using ISRestoreNonlocalIndices().
1446: The same scalability considerations as those for ISGetTotalIndices
1447: apply here.
1449: .seealso: ISGetTotalIndices(), ISRestoreNonlocalIndices(), ISGetSize(), ISGetLocalSize().
1450: @*/
1451: PetscErrorCode ISGetNonlocalIndices(IS is, const PetscInt *indices[])
1452: {
1454: PetscMPIInt size;
1455: PetscInt n, N;
1460: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
1461: if (size == 1) *indices = NULL;
1462: else {
1463: if (!is->total) {
1464: ISGatherTotal_Private(is);
1465: }
1466: ISGetLocalSize(is,&n);
1467: ISGetSize(is,&N);
1468: PetscMalloc1(N-n, &(is->nonlocal));
1469: PetscArraycpy(is->nonlocal, is->total, is->local_offset);
1470: PetscArraycpy(is->nonlocal+is->local_offset, is->total+is->local_offset+n,N - is->local_offset - n);
1471: *indices = is->nonlocal;
1472: }
1473: return(0);
1474: }
1476: /*@C
1477: ISRestoreTotalIndices - Restore the index array obtained with ISGetNonlocalIndices().
1479: Not Collective.
1481: Input Parameter:
1482: + is - the index set
1483: - indices - index array; must be the array obtained with ISGetNonlocalIndices()
1485: Level: intermediate
1487: .seealso: ISGetTotalIndices(), ISGetNonlocalIndices(), ISRestoreTotalIndices()
1488: @*/
1489: PetscErrorCode ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
1490: {
1494: 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.");
1495: return(0);
1496: }
1498: /*@
1499: ISGetNonlocalIS - Gather all nonlocal indices for this IS and present
1500: them as another sequential index set.
1502: Collective on IS
1504: Input Parameter:
1505: . is - the index set
1507: Output Parameter:
1508: . complement - sequential IS with indices identical to the result of
1509: ISGetNonlocalIndices()
1511: Level: intermediate
1513: Notes:
1514: complement represents the result of ISGetNonlocalIndices as an IS.
1515: Therefore scalability issues similar to ISGetNonlocalIndices apply.
1516: The resulting IS must be restored using ISRestoreNonlocalIS().
1518: .seealso: ISGetNonlocalIndices(), ISRestoreNonlocalIndices(), ISAllGather(), ISGetSize()
1519: @*/
1520: PetscErrorCode ISGetNonlocalIS(IS is, IS *complement)
1521: {
1527: /* Check if the complement exists already. */
1528: if (is->complement) {
1529: *complement = is->complement;
1530: PetscObjectReference((PetscObject)(is->complement));
1531: } else {
1532: PetscInt N, n;
1533: const PetscInt *idx;
1534: ISGetSize(is, &N);
1535: ISGetLocalSize(is,&n);
1536: ISGetNonlocalIndices(is, &idx);
1537: ISCreateGeneral(PETSC_COMM_SELF, N-n,idx, PETSC_USE_POINTER, &(is->complement));
1538: PetscObjectReference((PetscObject)is->complement);
1539: *complement = is->complement;
1540: }
1541: return(0);
1542: }
1544: /*@
1545: ISRestoreNonlocalIS - Restore the IS obtained with ISGetNonlocalIS().
1547: Not collective.
1549: Input Parameter:
1550: + is - the index set
1551: - complement - index set of is's nonlocal indices
1553: Level: intermediate
1555: .seealso: ISGetNonlocalIS(), ISGetNonlocalIndices(), ISRestoreNonlocalIndices()
1556: @*/
1557: PetscErrorCode ISRestoreNonlocalIS(IS is, IS *complement)
1558: {
1560: PetscInt refcnt;
1565: if (*complement != is->complement) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Complement IS being restored was not obtained with ISGetNonlocalIS()");
1566: PetscObjectGetReference((PetscObject)(is->complement), &refcnt);
1567: if (refcnt <= 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Duplicate call to ISRestoreNonlocalIS() detected");
1568: PetscObjectDereference((PetscObject)(is->complement));
1569: return(0);
1570: }
1572: /*@C
1573: ISViewFromOptions - View from Options
1575: Collective on IS
1577: Input Parameters:
1578: + A - the index set
1579: . obj - Optional object
1580: - name - command line option
1582: Level: intermediate
1583: .seealso: IS, ISView, PetscObjectViewFromOptions(), ISCreate()
1584: @*/
1585: PetscErrorCode ISViewFromOptions(IS A,PetscObject obj,const char name[])
1586: {
1591: PetscObjectViewFromOptions((PetscObject)A,obj,name);
1592: return(0);
1593: }
1595: /*@C
1596: ISView - Displays an index set.
1598: Collective on IS
1600: Input Parameters:
1601: + is - the index set
1602: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
1604: Level: intermediate
1606: .seealso: PetscViewerASCIIOpen()
1607: @*/
1608: PetscErrorCode ISView(IS is,PetscViewer viewer)
1609: {
1614: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is),&viewer);}
1618: PetscObjectPrintClassNamePrefixType((PetscObject)is,viewer);
1619: PetscLogEventBegin(IS_View,is,viewer,0,0);
1620: (*is->ops->view)(is,viewer);
1621: PetscLogEventEnd(IS_View,is,viewer,0,0);
1622: return(0);
1623: }
1625: /*@
1626: ISLoad - Loads a vector that has been stored in binary or HDF5 format with ISView().
1628: Collective on PetscViewer
1630: Input Parameters:
1631: + is - the newly loaded vector, this needs to have been created with ISCreate() or some related function before a call to ISLoad().
1632: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or HDF5 file viewer, obtained from PetscViewerHDF5Open()
1634: Level: intermediate
1636: Notes:
1637: IF using HDF5, you must assign the IS the same name as was used in the IS
1638: that was stored in the file using PetscObjectSetName(). Otherwise you will
1639: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"
1641: .seealso: PetscViewerBinaryOpen(), ISView(), MatLoad(), VecLoad()
1642: @*/
1643: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
1644: {
1645: PetscBool isbinary, ishdf5;
1652: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
1653: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
1654: if (!isbinary && !ishdf5) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1655: if (!((PetscObject)is)->type_name) {ISSetType(is, ISGENERAL);}
1656: PetscLogEventBegin(IS_Load,is,viewer,0,0);
1657: (*is->ops->load)(is, viewer);
1658: PetscLogEventEnd(IS_Load,is,viewer,0,0);
1659: return(0);
1660: }
1662: /*@
1663: ISSort - Sorts the indices of an index set.
1665: Collective on IS
1667: Input Parameters:
1668: . is - the index set
1670: Level: intermediate
1672: .seealso: ISSortRemoveDups(), ISSorted()
1673: @*/
1674: PetscErrorCode ISSort(IS is)
1675: {
1680: (*is->ops->sort)(is);
1681: ISSetInfo(is,IS_SORTED,IS_LOCAL,is->info_permanent[IS_LOCAL][IS_SORTED],PETSC_TRUE);
1682: return(0);
1683: }
1685: /*@
1686: ISSortRemoveDups - Sorts the indices of an index set, removing duplicates.
1688: Collective on IS
1690: Input Parameters:
1691: . is - the index set
1693: Level: intermediate
1695: .seealso: ISSort(), ISSorted()
1696: @*/
1697: PetscErrorCode ISSortRemoveDups(IS is)
1698: {
1703: ISClearInfoCache(is,PETSC_FALSE);
1704: (*is->ops->sortremovedups)(is);
1705: ISSetInfo(is,IS_SORTED,IS_LOCAL,is->info_permanent[IS_LOCAL][IS_SORTED],PETSC_TRUE);
1706: ISSetInfo(is,IS_UNIQUE,IS_LOCAL,is->info_permanent[IS_LOCAL][IS_UNIQUE],PETSC_TRUE);
1707: return(0);
1708: }
1710: /*@
1711: ISToGeneral - Converts an IS object of any type to ISGENERAL type
1713: Collective on IS
1715: Input Parameters:
1716: . is - the index set
1718: Level: intermediate
1720: .seealso: ISSorted()
1721: @*/
1722: PetscErrorCode ISToGeneral(IS is)
1723: {
1728: if (is->ops->togeneral) {
1729: (*is->ops->togeneral)(is);
1730: } else SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Not written for this type %s",((PetscObject)is)->type_name);
1731: return(0);
1732: }
1734: /*@
1735: ISSorted - Checks the indices to determine whether they have been sorted.
1737: Collective on IS
1739: Input Parameter:
1740: . is - the index set
1742: Output Parameter:
1743: . flg - output flag, either PETSC_TRUE if the index set is sorted,
1744: or PETSC_FALSE otherwise.
1746: Notes:
1747: For parallel IS objects this only indicates if the local part of the IS
1748: is sorted. So some processors may return PETSC_TRUE while others may
1749: return PETSC_FALSE.
1751: Level: intermediate
1753: .seealso: ISSort(), ISSortRemoveDups()
1754: @*/
1755: PetscErrorCode ISSorted(IS is,PetscBool *flg)
1756: {
1762: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,flg);
1763: return(0);
1764: }
1766: /*@
1767: ISDuplicate - Creates a duplicate copy of an index set.
1769: Collective on IS
1771: Input Parmeters:
1772: . is - the index set
1774: Output Parameters:
1775: . isnew - the copy of the index set
1777: Level: beginner
1779: .seealso: ISCreateGeneral(), ISCopy()
1780: @*/
1781: PetscErrorCode ISDuplicate(IS is,IS *newIS)
1782: {
1788: (*is->ops->duplicate)(is,newIS);
1789: ISCopyInfo(is,*newIS);
1790: return(0);
1791: }
1793: /*@
1794: ISCopy - Copies an index set.
1796: Collective on IS
1798: Input Parmeters:
1799: . is - the index set
1801: Output Parameters:
1802: . isy - the copy of the index set
1804: Level: beginner
1806: .seealso: ISDuplicate()
1807: @*/
1808: PetscErrorCode ISCopy(IS is,IS isy)
1809: {
1816: if (is == isy) return(0);
1817: ISCopyInfo(is,isy);
1818: isy->max = is->max;
1819: isy->min = is->min;
1820: (*is->ops->copy)(is,isy);
1821: return(0);
1822: }
1824: /*@
1825: ISOnComm - Split a parallel IS on subcomms (usually self) or concatenate index sets on subcomms into a parallel index set
1827: Collective on IS
1829: Input Arguments:
1830: + is - index set
1831: . comm - communicator for new index set
1832: - mode - copy semantics, PETSC_USE_POINTER for no-copy if possible, otherwise PETSC_COPY_VALUES
1834: Output Arguments:
1835: . newis - new IS on comm
1837: Level: advanced
1839: Notes:
1840: It is usually desirable to create a parallel IS and look at the local part when necessary.
1842: This function is useful if serial ISs must be created independently, or to view many
1843: logically independent serial ISs.
1845: The input IS must have the same type on every process.
1846: @*/
1847: PetscErrorCode ISOnComm(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
1848: {
1850: PetscMPIInt match;
1855: MPI_Comm_compare(PetscObjectComm((PetscObject)is),comm,&match);
1856: if (mode != PETSC_COPY_VALUES && (match == MPI_IDENT || match == MPI_CONGRUENT)) {
1857: PetscObjectReference((PetscObject)is);
1858: *newis = is;
1859: } else {
1860: (*is->ops->oncomm)(is,comm,mode,newis);
1861: }
1862: return(0);
1863: }
1865: /*@
1866: ISSetBlockSize - informs an index set that it has a given block size
1868: Logicall Collective on IS
1870: Input Arguments:
1871: + is - index set
1872: - bs - block size
1874: Level: intermediate
1876: Notes:
1877: This is much like the block size for Vecs. It indicates that one can think of the indices as
1878: being in a collection of equal size blocks. For ISBlock() these collections of blocks are all contiquous
1879: within a block but this is not the case for other IS.
1880: ISBlockGetIndices() only works for ISBlock IS, not others.
1882: .seealso: ISGetBlockSize(), ISCreateBlock(), ISBlockGetIndices(),
1883: @*/
1884: PetscErrorCode ISSetBlockSize(IS is,PetscInt bs)
1885: {
1891: if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_OUTOFRANGE,"Block size %D, must be positive",bs);
1892: (*is->ops->setblocksize)(is,bs);
1893: return(0);
1894: }
1896: /*@
1897: ISGetBlockSize - Returns the number of elements in a block.
1899: Not Collective
1901: Input Parameter:
1902: . is - the index set
1904: Output Parameter:
1905: . size - the number of elements in a block
1907: Level: intermediate
1909: Notes:
1910: This is much like the block size for Vecs. It indicates that one can think of the indices as
1911: being in a collection of equal size blocks. For ISBlock() these collections of blocks are all contiquous
1912: within a block but this is not the case for other IS.
1913: ISBlockGetIndices() only works for ISBlock IS, not others.
1915: .seealso: ISBlockGetSize(), ISGetSize(), ISCreateBlock(), ISSetBlockSize()
1916: @*/
1917: PetscErrorCode ISGetBlockSize(IS is,PetscInt *size)
1918: {
1922: PetscLayoutGetBlockSize(is->map, size);
1923: return(0);
1924: }
1926: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[])
1927: {
1929: PetscInt len,i;
1930: const PetscInt *ptr;
1933: ISGetLocalSize(is,&len);
1934: ISGetIndices(is,&ptr);
1935: for (i=0; i<len; i++) idx[i] = ptr[i];
1936: ISRestoreIndices(is,&ptr);
1937: return(0);
1938: }
1940: /*MC
1941: ISGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1942: The users should call ISRestoreIndicesF90() after having looked at the
1943: indices. The user should NOT change the indices.
1945: Synopsis:
1946: ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1948: Not collective
1950: Input Parameter:
1951: . x - index set
1953: Output Parameters:
1954: + xx_v - the Fortran90 pointer to the array
1955: - ierr - error code
1957: Example of Usage:
1958: .vb
1959: PetscInt, pointer xx_v(:)
1960: ....
1961: call ISGetIndicesF90(x,xx_v,ierr)
1962: a = xx_v(3)
1963: call ISRestoreIndicesF90(x,xx_v,ierr)
1964: .ve
1966: Level: intermediate
1968: .seealso: ISRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices()
1970: M*/
1972: /*MC
1973: ISRestoreIndicesF90 - Restores an index set to a usable state after
1974: a call to ISGetIndicesF90().
1976: Synopsis:
1977: ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1979: Not collective
1981: Input Parameters:
1982: + x - index set
1983: - xx_v - the Fortran90 pointer to the array
1985: Output Parameter:
1986: . ierr - error code
1988: Example of Usage:
1989: .vb
1990: PetscInt, pointer xx_v(:)
1991: ....
1992: call ISGetIndicesF90(x,xx_v,ierr)
1993: a = xx_v(3)
1994: call ISRestoreIndicesF90(x,xx_v,ierr)
1995: .ve
1997: Level: intermediate
1999: .seealso: ISGetIndicesF90(), ISGetIndices(), ISRestoreIndices()
2001: M*/
2003: /*MC
2004: ISBlockGetIndicesF90 - Accesses the elements of an index set from Fortran90.
2005: The users should call ISBlockRestoreIndicesF90() after having looked at the
2006: indices. The user should NOT change the indices.
2008: Synopsis:
2009: ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
2011: Not collective
2013: Input Parameter:
2014: . x - index set
2016: Output Parameters:
2017: + xx_v - the Fortran90 pointer to the array
2018: - ierr - error code
2019: Example of Usage:
2020: .vb
2021: PetscInt, pointer xx_v(:)
2022: ....
2023: call ISBlockGetIndicesF90(x,xx_v,ierr)
2024: a = xx_v(3)
2025: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
2026: .ve
2028: Level: intermediate
2030: .seealso: ISBlockRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices(),
2031: ISRestoreIndices()
2033: M*/
2035: /*MC
2036: ISBlockRestoreIndicesF90 - Restores an index set to a usable state after
2037: a call to ISBlockGetIndicesF90().
2039: Synopsis:
2040: ISBlockRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
2042: Not Collective
2044: Input Parameters:
2045: + x - index set
2046: - xx_v - the Fortran90 pointer to the array
2048: Output Parameter:
2049: . ierr - error code
2051: Example of Usage:
2052: .vb
2053: PetscInt, pointer xx_v(:)
2054: ....
2055: call ISBlockGetIndicesF90(x,xx_v,ierr)
2056: a = xx_v(3)
2057: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
2058: .ve
2060: Notes:
2061: Not yet supported for all F90 compilers
2063: Level: intermediate
2065: .seealso: ISBlockGetIndicesF90(), ISGetIndices(), ISRestoreIndices(), ISRestoreIndicesF90()
2067: M*/