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 the non-negative entries of an index set in a contiguous way, starting from 0.
17: Collective on IS
19: Input Parameters:
20: + subset - the index set
21: - subset_mult - the multiplicity 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: Notes: All negative entries are mapped to -1. Indices with non positive multiplicities are skipped.
29: Level: intermediate
31: .seealso:
32: @*/
33: PetscErrorCode ISRenumber(IS subset, IS subset_mult, PetscInt *N, IS *subset_n)
34: {
35: PetscSF sf;
36: PetscLayout map;
37: const PetscInt *idxs, *idxs_mult = NULL;
38: PetscInt *leaf_data,*root_data,*gidxs,*ilocal,*ilocalneg;
39: PetscInt N_n,n,i,lbounds[2],gbounds[2],Nl,ibs;
40: PetscInt n_n,nlocals,start,first_index,npos,nneg;
41: PetscMPIInt commsize;
42: PetscBool first_found,isblock;
47: else if (!subset_n) return 0;
48: ISGetLocalSize(subset,&n);
49: if (subset_mult) {
50: ISGetLocalSize(subset_mult,&i);
52: }
53: /* create workspace layout for computing global indices of subset */
54: PetscMalloc1(n,&ilocal);
55: PetscMalloc1(n,&ilocalneg);
56: ISGetIndices(subset,&idxs);
57: ISGetBlockSize(subset,&ibs);
58: PetscObjectTypeCompare((PetscObject)subset,ISBLOCK,&isblock);
59: if (subset_mult) ISGetIndices(subset_mult,&idxs_mult);
60: lbounds[0] = PETSC_MAX_INT;
61: lbounds[1] = PETSC_MIN_INT;
62: for (i=0,npos=0,nneg=0;i<n;i++) {
63: if (idxs[i] < 0) { ilocalneg[nneg++] = i; continue; }
64: if (idxs[i] < lbounds[0]) lbounds[0] = idxs[i];
65: if (idxs[i] > lbounds[1]) lbounds[1] = idxs[i];
66: ilocal[npos++] = i;
67: }
68: if (npos == n) {
69: PetscFree(ilocal);
70: PetscFree(ilocalneg);
71: }
73: /* create sf : leaf_data == multiplicity of indexes, root data == global index in layout */
74: PetscMalloc1(n,&leaf_data);
75: for (i=0;i<n;i++) leaf_data[i] = idxs_mult ? PetscMax(idxs_mult[i],0) : 1;
77: /* local size of new subset */
78: n_n = 0;
79: for (i=0;i<n;i++) n_n += leaf_data[i];
80: if (ilocalneg) for (i=0;i<nneg;i++) leaf_data[ilocalneg[i]] = 0;
81: PetscFree(ilocalneg);
82: PetscMalloc1(PetscMax(n_n,n),&gidxs); /* allocating extra space to reuse gidxs */
83: /* check for early termination (all negative) */
84: PetscGlobalMinMaxInt(PetscObjectComm((PetscObject)subset),lbounds,gbounds);
85: if (gbounds[1] < gbounds[0]) {
86: if (N) *N = 0;
87: if (subset_n) { /* all negative */
88: for (i=0;i<n_n;i++) gidxs[i] = -1;
89: ISCreateGeneral(PetscObjectComm((PetscObject)subset),n_n,gidxs,PETSC_COPY_VALUES,subset_n);
90: }
91: PetscFree(leaf_data);
92: PetscFree(gidxs);
93: ISRestoreIndices(subset,&idxs);
94: if (subset_mult) ISRestoreIndices(subset_mult,&idxs_mult);
95: PetscFree(ilocal);
96: PetscFree(ilocalneg);
97: return 0;
98: }
100: /* split work */
101: N_n = gbounds[1] - gbounds[0] + 1;
102: PetscLayoutCreate(PetscObjectComm((PetscObject)subset),&map);
103: PetscLayoutSetBlockSize(map,1);
104: PetscLayoutSetSize(map,N_n);
105: PetscLayoutSetUp(map);
106: PetscLayoutGetLocalSize(map,&Nl);
108: /* global indexes in layout */
109: for (i=0;i<npos;i++) gidxs[i] = (ilocal ? idxs[ilocal[i]] : idxs[i]) - gbounds[0];
110: ISRestoreIndices(subset,&idxs);
111: PetscSFCreate(PetscObjectComm((PetscObject)subset),&sf);
112: PetscSFSetGraphLayout(sf,map,npos,ilocal,PETSC_USE_POINTER,gidxs);
113: PetscLayoutDestroy(&map);
115: /* reduce from leaves to roots */
116: PetscCalloc1(Nl,&root_data);
117: PetscSFReduceBegin(sf,MPIU_INT,leaf_data,root_data,MPI_MAX);
118: PetscSFReduceEnd(sf,MPIU_INT,leaf_data,root_data,MPI_MAX);
120: /* count indexes in local part of layout */
121: nlocals = 0;
122: first_index = -1;
123: first_found = PETSC_FALSE;
124: for (i=0;i<Nl;i++) {
125: if (!first_found && root_data[i]) {
126: first_found = PETSC_TRUE;
127: first_index = i;
128: }
129: nlocals += root_data[i];
130: }
132: /* cumulative of number of indexes and size of subset without holes */
133: #if defined(PETSC_HAVE_MPI_EXSCAN)
134: start = 0;
135: MPI_Exscan(&nlocals,&start,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)subset));
136: #else
137: MPI_Scan(&nlocals,&start,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)subset));
138: start = start-nlocals;
139: #endif
141: if (N) { /* compute total size of new subset if requested */
142: *N = start + nlocals;
143: MPI_Comm_size(PetscObjectComm((PetscObject)subset),&commsize);
144: MPI_Bcast(N,1,MPIU_INT,commsize-1,PetscObjectComm((PetscObject)subset));
145: }
147: if (!subset_n) {
148: PetscFree(gidxs);
149: PetscSFDestroy(&sf);
150: PetscFree(leaf_data);
151: PetscFree(root_data);
152: PetscFree(ilocal);
153: if (subset_mult) ISRestoreIndices(subset_mult,&idxs_mult);
154: return 0;
155: }
157: /* adapt root data with cumulative */
158: if (first_found) {
159: PetscInt old_index;
161: root_data[first_index] += start;
162: old_index = first_index;
163: for (i=first_index+1;i<Nl;i++) {
164: if (root_data[i]) {
165: root_data[i] += root_data[old_index];
166: old_index = i;
167: }
168: }
169: }
171: /* from roots to leaves */
172: PetscSFBcastBegin(sf,MPIU_INT,root_data,leaf_data,MPI_REPLACE);
173: PetscSFBcastEnd(sf,MPIU_INT,root_data,leaf_data,MPI_REPLACE);
174: PetscSFDestroy(&sf);
176: /* create new IS with global indexes without holes */
177: for (i=0;i<n_n;i++) gidxs[i] = -1;
178: if (subset_mult) {
179: PetscInt cum;
181: isblock = PETSC_FALSE;
182: for (i=0,cum=0;i<n;i++) for (PetscInt j=0;j<idxs_mult[i];j++) gidxs[cum++] = leaf_data[i] - idxs_mult[i] + j;
183: } else for (i=0;i<n;i++) gidxs[i] = leaf_data[i]-1;
185: if (isblock) {
186: if (ibs > 1) for (i=0;i<n_n/ibs;i++) gidxs[i] = gidxs[i*ibs]/ibs;
187: ISCreateBlock(PetscObjectComm((PetscObject)subset),ibs,n_n/ibs,gidxs,PETSC_COPY_VALUES,subset_n);
188: } else {
189: ISCreateGeneral(PetscObjectComm((PetscObject)subset),n_n,gidxs,PETSC_COPY_VALUES,subset_n);
190: }
191: if (subset_mult) ISRestoreIndices(subset_mult,&idxs_mult);
192: PetscFree(gidxs);
193: PetscFree(leaf_data);
194: PetscFree(root_data);
195: PetscFree(ilocal);
196: return 0;
197: }
199: /*@
200: ISCreateSubIS - Create a sub index set from a global index set selecting some components.
202: Collective on IS
204: Input Parameters:
205: + is - the index set
206: - comps - which components we will extract from is
208: Output Parameters:
209: . subis - the new sub index set
211: Level: intermediate
213: Example usage:
214: We have an index set (is) living on 3 processes with the following values:
215: | 4 9 0 | 2 6 7 | 10 11 1|
216: and another index set (comps) used to indicate which components of is we want to take,
217: | 7 5 | 1 2 | 0 4|
218: The output index set (subis) should look like:
219: | 11 7 | 9 0 | 4 6|
221: .seealso: VecGetSubVector(), MatCreateSubMatrix()
222: @*/
223: PetscErrorCode ISCreateSubIS(IS is,IS comps,IS *subis)
224: {
225: PetscSF sf;
226: const PetscInt *is_indices,*comps_indices;
227: PetscInt *subis_indices,nroots,nleaves,*mine,i,lidx;
228: PetscMPIInt owner;
229: PetscSFNode *remote;
230: MPI_Comm comm;
236: PetscObjectGetComm((PetscObject)is, &comm);
237: ISGetLocalSize(comps,&nleaves);
238: ISGetLocalSize(is,&nroots);
239: PetscMalloc1(nleaves,&remote);
240: PetscMalloc1(nleaves,&mine);
241: ISGetIndices(comps,&comps_indices);
242: /*
243: * Construct a PetscSF in which "is" data serves as roots and "subis" is leaves.
244: * Root data are sent to leaves using PetscSFBcast().
245: * */
246: for (i=0; i<nleaves; i++) {
247: mine[i] = i;
248: /* Connect a remote root with the current leaf. The value on the remote root
249: * will be received by the current local leaf.
250: * */
251: owner = -1;
252: lidx = -1;
253: PetscLayoutFindOwnerIndex(is->map,comps_indices[i],&owner,&lidx);
254: remote[i].rank = owner;
255: remote[i].index = lidx;
256: }
257: ISRestoreIndices(comps,&comps_indices);
258: PetscSFCreate(comm,&sf);
259: PetscSFSetFromOptions(sf);\
260: PetscSFSetGraph(sf,nroots,nleaves,mine,PETSC_OWN_POINTER,remote,PETSC_OWN_POINTER);
262: PetscMalloc1(nleaves,&subis_indices);
263: ISGetIndices(is, &is_indices);
264: PetscSFBcastBegin(sf,MPIU_INT,is_indices,subis_indices,MPI_REPLACE);
265: PetscSFBcastEnd(sf,MPIU_INT,is_indices,subis_indices,MPI_REPLACE);
266: ISRestoreIndices(is,&is_indices);
267: PetscSFDestroy(&sf);
268: ISCreateGeneral(comm,nleaves,subis_indices,PETSC_OWN_POINTER,subis);
269: return 0;
270: }
272: /*@
273: ISClearInfoCache - clear the cache of computed index set properties
275: Not collective
277: Input Parameters:
278: + is - the index set
279: - clear_permanent_local - whether to remove the permanent status of local properties
281: NOTE: because all processes must agree on the global permanent status of a property,
282: the permanent status can only be changed with ISSetInfo(), because this routine is not collective
284: Level: developer
286: .seealso: ISInfo, ISInfoType, ISSetInfo(), ISClearInfoCache()
288: @*/
289: PetscErrorCode ISClearInfoCache(IS is, PetscBool clear_permanent_local)
290: {
291: PetscInt i, j;
295: for (i = 0; i < IS_INFO_MAX; i++) {
296: if (clear_permanent_local) is->info_permanent[IS_LOCAL][i] = PETSC_FALSE;
297: for (j = 0; j < 2; j++) {
298: if (!is->info_permanent[j][i]) is->info[j][i] = IS_INFO_UNKNOWN;
299: }
300: }
301: return 0;
302: }
304: static PetscErrorCode ISSetInfo_Internal(IS is, ISInfo info, ISInfoType type, ISInfoBool ipermanent, PetscBool flg)
305: {
306: ISInfoBool iflg = flg ? IS_INFO_TRUE : IS_INFO_FALSE;
307: PetscInt itype = (type == IS_LOCAL) ? 0 : 1;
308: PetscBool permanent_set = (ipermanent == IS_INFO_UNKNOWN) ? PETSC_FALSE : PETSC_TRUE;
309: PetscBool permanent = (ipermanent == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
311: /* set this property */
312: is->info[itype][(int)info] = iflg;
313: if (permanent_set) is->info_permanent[itype][(int)info] = permanent;
314: /* set implications */
315: switch (info) {
316: case IS_SORTED:
317: if (flg && type == IS_GLOBAL) { /* an array that is globally sorted is also locally sorted */
318: is->info[IS_LOCAL][(int)info] = IS_INFO_TRUE;
319: /* global permanence implies local permanence */
320: if (permanent_set && permanent) is->info_permanent[IS_LOCAL][(int)info] = PETSC_TRUE;
321: }
322: if (!flg) { /* if an array is not sorted, it cannot be an interval or the identity */
323: is->info[itype][IS_INTERVAL] = IS_INFO_FALSE;
324: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
325: if (permanent_set) {
326: is->info_permanent[itype][IS_INTERVAL] = permanent;
327: is->info_permanent[itype][IS_IDENTITY] = permanent;
328: }
329: }
330: break;
331: case IS_UNIQUE:
332: if (flg && type == IS_GLOBAL) { /* an array that is globally unique is also locally unique */
333: is->info[IS_LOCAL][(int)info] = IS_INFO_TRUE;
334: /* global permanence implies local permanence */
335: if (permanent_set && permanent) is->info_permanent[IS_LOCAL][(int)info] = PETSC_TRUE;
336: }
337: if (!flg) { /* if an array is not unique, it cannot be a permutation, and interval, or the identity */
338: is->info[itype][IS_PERMUTATION] = IS_INFO_FALSE;
339: is->info[itype][IS_INTERVAL] = IS_INFO_FALSE;
340: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
341: if (permanent_set) {
342: is->info_permanent[itype][IS_PERMUTATION] = permanent;
343: is->info_permanent[itype][IS_INTERVAL] = permanent;
344: is->info_permanent[itype][IS_IDENTITY] = permanent;
345: }
346: }
347: break;
348: case IS_PERMUTATION:
349: if (flg) { /* an array that is a permutation is unique and is unique locally */
350: is->info[itype][IS_UNIQUE] = IS_INFO_TRUE;
351: is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
352: if (permanent_set && permanent) {
353: is->info_permanent[itype][IS_UNIQUE] = PETSC_TRUE;
354: is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
355: }
356: } else { /* an array that is not a permutation cannot be the identity */
357: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
358: if (permanent_set) is->info_permanent[itype][IS_IDENTITY] = permanent;
359: }
360: break;
361: case IS_INTERVAL:
362: if (flg) { /* an array that is an interval is sorted and unique */
363: is->info[itype][IS_SORTED] = IS_INFO_TRUE;
364: is->info[IS_LOCAL][IS_SORTED] = IS_INFO_TRUE;
365: is->info[itype][IS_UNIQUE] = IS_INFO_TRUE;
366: is->info[IS_LOCAL][IS_UNIQUE] = 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: }
373: } else { /* an array that is not an interval cannot be the identity */
374: is->info[itype][IS_IDENTITY] = IS_INFO_FALSE;
375: if (permanent_set) is->info_permanent[itype][IS_IDENTITY] = permanent;
376: }
377: break;
378: case IS_IDENTITY:
379: if (flg) { /* an array that is the identity is sorted, unique, an interval, and a permutation */
380: is->info[itype][IS_SORTED] = IS_INFO_TRUE;
381: is->info[IS_LOCAL][IS_SORTED] = IS_INFO_TRUE;
382: is->info[itype][IS_UNIQUE] = IS_INFO_TRUE;
383: is->info[IS_LOCAL][IS_UNIQUE] = IS_INFO_TRUE;
384: is->info[itype][IS_PERMUTATION] = IS_INFO_TRUE;
385: is->info[itype][IS_INTERVAL] = IS_INFO_TRUE;
386: is->info[IS_LOCAL][IS_INTERVAL] = IS_INFO_TRUE;
387: if (permanent_set && permanent) {
388: is->info_permanent[itype][IS_SORTED] = PETSC_TRUE;
389: is->info_permanent[IS_LOCAL][IS_SORTED] = PETSC_TRUE;
390: is->info_permanent[itype][IS_UNIQUE] = PETSC_TRUE;
391: is->info_permanent[IS_LOCAL][IS_UNIQUE] = PETSC_TRUE;
392: is->info_permanent[itype][IS_PERMUTATION] = PETSC_TRUE;
393: is->info_permanent[itype][IS_INTERVAL] = PETSC_TRUE;
394: is->info_permanent[IS_LOCAL][IS_INTERVAL] = PETSC_TRUE;
395: }
396: }
397: break;
398: default:
400: else SETERRQ(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
401: }
402: return 0;
403: }
405: /*@
406: ISSetInfo - Set known information about an index set.
408: Logically Collective on IS if type is IS_GLOBAL
410: Input Parameters:
411: + is - the index set
412: . info - describing a property of the index set, one of those listed below,
413: . type - IS_LOCAL if the information describes the local portion of the index set,
414: IS_GLOBAL if it describes the whole index set
415: . permanent - PETSC_TRUE if it is known that the property will persist through changes to the index set, PETSC_FALSE otherwise
416: If the user sets a property as permanently known, it will bypass computation of that property
417: - flg - set the described property as true (PETSC_TRUE) or false (PETSC_FALSE)
419: Info Describing IS Structure:
420: + IS_SORTED - the [local part of the] index set is sorted in ascending order
421: . IS_UNIQUE - each entry in the [local part of the] index set is unique
422: . 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
423: . IS_INTERVAL - the [local part of the] index set is equal to a contiguous range of integers {f, f + 1, ..., f + N-1}
424: - IS_IDENTITY - the [local part of the] index set is equal to the integers {0, 1, ..., N-1}
426: Notes:
427: If type is IS_GLOBAL, all processes that share the index set must pass the same value in flg
429: It is possible to set a property with ISSetInfo() that contradicts what would be previously computed with ISGetInfo()
431: Level: advanced
433: .seealso: ISInfo, ISInfoType, IS
435: @*/
436: PetscErrorCode ISSetInfo(IS is, ISInfo info, ISInfoType type, PetscBool permanent, PetscBool flg)
437: {
438: MPI_Comm comm, errcomm;
439: PetscMPIInt size;
443: comm = PetscObjectComm((PetscObject)is);
444: if (type == IS_GLOBAL) {
448: errcomm = comm;
449: } else {
450: errcomm = PETSC_COMM_SELF;
451: }
455: MPI_Comm_size(comm, &size);
456: /* do not use global values if size == 1: it makes it easier to keep the implications straight */
457: if (size == 1) type = IS_LOCAL;
458: ISSetInfo_Internal(is, info, type, permanent ? IS_INFO_TRUE : IS_INFO_FALSE, flg);
459: return 0;
460: }
462: static PetscErrorCode ISGetInfo_Sorted(IS is, ISInfoType type, PetscBool *flg)
463: {
464: MPI_Comm comm;
465: PetscMPIInt size, rank;
467: comm = PetscObjectComm((PetscObject)is);
468: MPI_Comm_size(comm, &size);
469: MPI_Comm_size(comm, &rank);
470: if (type == IS_GLOBAL && is->ops->sortedglobal) {
471: (*is->ops->sortedglobal)(is,flg);
472: } else {
473: PetscBool sortedLocal = PETSC_FALSE;
475: /* determine if the array is locally sorted */
476: if (type == IS_GLOBAL && size > 1) {
477: /* call ISGetInfo so that a cached value will be used if possible */
478: ISGetInfo(is, IS_SORTED, IS_LOCAL, PETSC_TRUE, &sortedLocal);
479: } else if (is->ops->sortedlocal) {
480: (*is->ops->sortedlocal)(is,&sortedLocal);
481: } else {
482: /* default: get the local indices and directly check */
483: const PetscInt *idx;
484: PetscInt n;
486: ISGetIndices(is, &idx);
487: ISGetLocalSize(is, &n);
488: PetscSortedInt(n, idx, &sortedLocal);
489: ISRestoreIndices(is, &idx);
490: }
492: if (type == IS_LOCAL || size == 1) {
493: *flg = sortedLocal;
494: } else {
495: MPI_Allreduce(&sortedLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
496: if (*flg) {
497: PetscInt n, min = PETSC_MAX_INT, max = PETSC_MIN_INT;
498: PetscInt maxprev;
500: ISGetLocalSize(is, &n);
501: if (n) ISGetMinMax(is, &min, &max);
502: maxprev = PETSC_MIN_INT;
503: MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm);
504: if (rank && (maxprev > min)) sortedLocal = PETSC_FALSE;
505: MPI_Allreduce(&sortedLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
506: }
507: }
508: }
509: return 0;
510: }
512: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[]);
514: static PetscErrorCode ISGetInfo_Unique(IS is, ISInfoType type, PetscBool *flg)
515: {
516: MPI_Comm comm;
517: PetscMPIInt size, rank;
518: PetscInt i;
520: comm = PetscObjectComm((PetscObject)is);
521: MPI_Comm_size(comm, &size);
522: MPI_Comm_size(comm, &rank);
523: if (type == IS_GLOBAL && is->ops->uniqueglobal) {
524: (*is->ops->uniqueglobal)(is,flg);
525: } else {
526: PetscBool uniqueLocal;
527: PetscInt n = -1;
528: PetscInt *idx = NULL;
530: /* determine if the array is locally unique */
531: if (type == IS_GLOBAL && size > 1) {
532: /* call ISGetInfo so that a cached value will be used if possible */
533: ISGetInfo(is, IS_UNIQUE, IS_LOCAL, PETSC_TRUE, &uniqueLocal);
534: } else if (is->ops->uniquelocal) {
535: (*is->ops->uniquelocal)(is,&uniqueLocal);
536: } else {
537: /* default: get the local indices and directly check */
538: uniqueLocal = PETSC_TRUE;
539: ISGetLocalSize(is, &n);
540: PetscMalloc1(n, &idx);
541: ISGetIndicesCopy(is, idx);
542: PetscIntSortSemiOrdered(n, idx);
543: for (i = 1; i < n; i++) if (idx[i] == idx[i-1]) break;
544: if (i < n) uniqueLocal = PETSC_FALSE;
545: }
547: PetscFree(idx);
548: if (type == IS_LOCAL || size == 1) {
549: *flg = uniqueLocal;
550: } else {
551: MPI_Allreduce(&uniqueLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
552: if (*flg) {
553: PetscInt min = PETSC_MAX_INT, max = PETSC_MIN_INT, maxprev;
555: if (!idx) {
556: ISGetLocalSize(is, &n);
557: PetscMalloc1(n, &idx);
558: ISGetIndicesCopy(is, idx);
559: }
560: PetscParallelSortInt(is->map, is->map, idx, idx);
561: if (n) {
562: min = idx[0];
563: max = idx[n - 1];
564: }
565: for (i = 1; i < n; i++) if (idx[i] == idx[i-1]) break;
566: if (i < n) uniqueLocal = PETSC_FALSE;
567: maxprev = PETSC_MIN_INT;
568: MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm);
569: if (rank && (maxprev == min)) uniqueLocal = PETSC_FALSE;
570: MPI_Allreduce(&uniqueLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
571: }
572: }
573: PetscFree(idx);
574: }
575: return 0;
576: }
578: static PetscErrorCode ISGetInfo_Permutation(IS is, ISInfoType type, PetscBool *flg)
579: {
580: MPI_Comm comm;
581: PetscMPIInt size, rank;
583: comm = PetscObjectComm((PetscObject)is);
584: MPI_Comm_size(comm, &size);
585: MPI_Comm_size(comm, &rank);
586: if (type == IS_GLOBAL && is->ops->permglobal) {
587: (*is->ops->permglobal)(is,flg);
588: } else if (type == IS_LOCAL && is->ops->permlocal) {
589: (*is->ops->permlocal)(is,flg);
590: } else {
591: PetscBool permLocal;
592: PetscInt n, i, rStart;
593: PetscInt *idx;
595: ISGetLocalSize(is, &n);
596: PetscMalloc1(n, &idx);
597: ISGetIndicesCopy(is, idx);
598: if (type == IS_GLOBAL) {
599: PetscParallelSortInt(is->map, is->map, idx, idx);
600: PetscLayoutGetRange(is->map, &rStart, NULL);
601: } else {
602: PetscIntSortSemiOrdered(n, idx);
603: rStart = 0;
604: }
605: permLocal = PETSC_TRUE;
606: for (i = 0; i < n; i++) {
607: if (idx[i] != rStart + i) break;
608: }
609: if (i < n) permLocal = PETSC_FALSE;
610: if (type == IS_LOCAL || size == 1) {
611: *flg = permLocal;
612: } else {
613: MPI_Allreduce(&permLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
614: }
615: PetscFree(idx);
616: }
617: return 0;
618: }
620: static PetscErrorCode ISGetInfo_Interval(IS is, ISInfoType type, PetscBool *flg)
621: {
622: MPI_Comm comm;
623: PetscMPIInt size, rank;
624: PetscInt i;
626: comm = PetscObjectComm((PetscObject)is);
627: MPI_Comm_size(comm, &size);
628: MPI_Comm_size(comm, &rank);
629: if (type == IS_GLOBAL && is->ops->intervalglobal) {
630: (*is->ops->intervalglobal)(is,flg);
631: } else {
632: PetscBool intervalLocal;
634: /* determine if the array is locally an interval */
635: if (type == IS_GLOBAL && size > 1) {
636: /* call ISGetInfo so that a cached value will be used if possible */
637: ISGetInfo(is, IS_INTERVAL, IS_LOCAL, PETSC_TRUE, &intervalLocal);
638: } else if (is->ops->intervallocal) {
639: (*is->ops->intervallocal)(is,&intervalLocal);
640: } else {
641: PetscInt n;
642: const PetscInt *idx;
643: /* default: get the local indices and directly check */
644: intervalLocal = PETSC_TRUE;
645: ISGetLocalSize(is, &n);
646: PetscMalloc1(n, &idx);
647: ISGetIndices(is, &idx);
648: for (i = 1; i < n; i++) if (idx[i] != idx[i-1] + 1) break;
649: if (i < n) intervalLocal = PETSC_FALSE;
650: ISRestoreIndices(is, &idx);
651: }
653: if (type == IS_LOCAL || size == 1) {
654: *flg = intervalLocal;
655: } else {
656: MPI_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
657: if (*flg) {
658: PetscInt n, min = PETSC_MAX_INT, max = PETSC_MIN_INT;
659: PetscInt maxprev;
661: ISGetLocalSize(is, &n);
662: if (n) ISGetMinMax(is, &min, &max);
663: maxprev = PETSC_MIN_INT;
664: MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm);
665: if (rank && n && (maxprev != min - 1)) intervalLocal = PETSC_FALSE;
666: MPI_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
667: }
668: }
669: }
670: return 0;
671: }
673: static PetscErrorCode ISGetInfo_Identity(IS is, ISInfoType type, PetscBool *flg)
674: {
675: MPI_Comm comm;
676: PetscMPIInt size, rank;
678: comm = PetscObjectComm((PetscObject)is);
679: MPI_Comm_size(comm, &size);
680: MPI_Comm_size(comm, &rank);
681: if (type == IS_GLOBAL && is->ops->intervalglobal) {
682: PetscBool isinterval;
684: (*is->ops->intervalglobal)(is,&isinterval);
685: *flg = PETSC_FALSE;
686: if (isinterval) {
687: PetscInt min;
689: ISGetMinMax(is, &min, NULL);
690: MPI_Bcast(&min, 1, MPIU_INT, 0, comm);
691: if (min == 0) *flg = PETSC_TRUE;
692: }
693: } else if (type == IS_LOCAL && is->ops->intervallocal) {
694: PetscBool isinterval;
696: (*is->ops->intervallocal)(is,&isinterval);
697: *flg = PETSC_FALSE;
698: if (isinterval) {
699: PetscInt min;
701: ISGetMinMax(is, &min, NULL);
702: if (min == 0) *flg = PETSC_TRUE;
703: }
704: } else {
705: PetscBool identLocal;
706: PetscInt n, i, rStart;
707: const PetscInt *idx;
709: ISGetLocalSize(is, &n);
710: ISGetIndices(is, &idx);
711: PetscLayoutGetRange(is->map, &rStart, NULL);
712: identLocal = PETSC_TRUE;
713: for (i = 0; i < n; i++) {
714: if (idx[i] != rStart + i) break;
715: }
716: if (i < n) identLocal = PETSC_FALSE;
717: if (type == IS_LOCAL || size == 1) {
718: *flg = identLocal;
719: } else {
720: MPI_Allreduce(&identLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
721: }
722: ISRestoreIndices(is, &idx);
723: }
724: return 0;
725: }
727: /*@
728: ISGetInfo - Determine whether an index set satisfies a given property
730: 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())
732: Input Parameters:
733: + is - the index set
734: . info - describing a property of the index set, one of those listed in the documentation of ISSetInfo()
735: . 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
736: - type - whether the property is local (IS_LOCAL) or global (IS_GLOBAL)
738: Output Parameter:
739: . flg - wheter the property is true (PETSC_TRUE) or false (PETSC_FALSE)
741: 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().
743: Level: advanced
745: .seealso: ISInfo, ISInfoType, ISSetInfo(), ISClearInfoCache()
747: @*/
748: PetscErrorCode ISGetInfo(IS is, ISInfo info, ISInfoType type, PetscBool compute, PetscBool *flg)
749: {
750: MPI_Comm comm, errcomm;
751: PetscMPIInt rank, size;
752: PetscInt itype;
753: PetscBool hasprop;
754: PetscBool infer;
758: comm = PetscObjectComm((PetscObject)is);
759: if (type == IS_GLOBAL) {
761: errcomm = comm;
762: } else {
763: errcomm = PETSC_COMM_SELF;
764: }
766: MPI_Comm_size(comm, &size);
767: MPI_Comm_rank(comm, &rank);
770: if (size == 1) type = IS_LOCAL;
771: itype = (type == IS_LOCAL) ? 0 : 1;
772: hasprop = PETSC_FALSE;
773: infer = PETSC_FALSE;
774: if (is->info_permanent[itype][(int)info]) {
775: hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
776: infer = PETSC_TRUE;
777: } else if ((itype == IS_LOCAL) && (is->info[IS_LOCAL][info] != IS_INFO_UNKNOWN)) {
778: /* we can cache local properties as long as we clear them when the IS changes */
779: /* NOTE: we only cache local values because there is no ISAssemblyBegin()/ISAssemblyEnd(),
780: so we have no way of knowing when a cached value has been invalidated by changes on a different process */
781: hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
782: infer = PETSC_TRUE;
783: } else if (compute) {
784: switch (info) {
785: case IS_SORTED:
786: ISGetInfo_Sorted(is, type, &hasprop);
787: break;
788: case IS_UNIQUE:
789: ISGetInfo_Unique(is, type, &hasprop);
790: break;
791: case IS_PERMUTATION:
792: ISGetInfo_Permutation(is, type, &hasprop);
793: break;
794: case IS_INTERVAL:
795: ISGetInfo_Interval(is, type, &hasprop);
796: break;
797: case IS_IDENTITY:
798: ISGetInfo_Identity(is, type, &hasprop);
799: break;
800: default:
801: SETERRQ(errcomm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
802: }
803: infer = PETSC_TRUE;
804: }
805: /* call ISSetInfo_Internal to keep all of the implications straight */
806: if (infer) ISSetInfo_Internal(is, info, type, IS_INFO_UNKNOWN, hasprop);
807: *flg = hasprop;
808: return 0;
809: }
811: static PetscErrorCode ISCopyInfo(IS source, IS dest)
812: {
813: PetscArraycpy(&dest->info[0], &source->info[0], 2);
814: PetscArraycpy(&dest->info_permanent[0], &source->info_permanent[0], 2);
815: return 0;
816: }
818: /*@
819: ISIdentity - Determines whether index set is the identity mapping.
821: Collective on IS
823: Input Parameters:
824: . is - the index set
826: Output Parameters:
827: . ident - PETSC_TRUE if an identity, else PETSC_FALSE
829: Level: intermediate
831: Note: If ISSetIdentity() (or ISSetInfo() for a permanent property) has been called,
832: ISIdentity() will return its answer without communication between processes, but
833: otherwise the output ident will be computed from ISGetInfo(),
834: which may require synchronization on the communicator of IS. To avoid this computation,
835: call ISGetInfo() directly with the compute flag set to PETSC_FALSE, and ident will be assumed false.
837: .seealso: ISSetIdentity(), ISGetInfo()
838: @*/
839: PetscErrorCode ISIdentity(IS is,PetscBool *ident)
840: {
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 Parameter:
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: {
866: ISSetInfo(is,IS_IDENTITY,IS_GLOBAL,PETSC_TRUE,PETSC_TRUE);
867: return 0;
868: }
870: /*@
871: ISContiguousLocal - Locates an index set with contiguous range within a global range, if possible
873: Not Collective
875: Input Parameters:
876: + is - the index set
877: . gstart - global start
878: - gend - global end
880: Output Parameters:
881: + start - start of contiguous block, as an offset from gstart
882: - contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
884: Level: developer
886: .seealso: ISGetLocalSize(), VecGetOwnershipRange()
887: @*/
888: PetscErrorCode ISContiguousLocal(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
889: {
893: *start = -1;
894: *contig = PETSC_FALSE;
895: if (is->ops->contiguous) {
896: (*is->ops->contiguous)(is,gstart,gend,start,contig);
897: }
898: return 0;
899: }
901: /*@
902: ISPermutation - PETSC_TRUE or PETSC_FALSE depending on whether the
903: index set has been declared to be a permutation.
905: Logically Collective on IS
907: Input Parameter:
908: . is - the index set
910: Output Parameter:
911: . perm - PETSC_TRUE if a permutation, else PETSC_FALSE
913: Level: intermediate
915: Note: If it is not alread known that the IS is a permutation (if ISSetPermutation()
916: or ISSetInfo() has not been called), this routine will not attempt to compute
917: whether the index set is a permutation and will assume perm is PETSC_FALSE.
918: To compute the value when it is not already known, use ISGetInfo() with
919: the compute flag set to PETSC_TRUE.
921: .seealso: ISSetPermutation(), ISGetInfo()
922: @*/
923: PetscErrorCode ISPermutation(IS is,PetscBool *perm)
924: {
927: ISGetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_FALSE,perm);
928: return 0;
929: }
931: /*@
932: ISSetPermutation - Informs the index set that it is a permutation.
934: Logically Collective on IS
936: Input Parameter:
937: . is - the index set
939: Level: intermediate
941: The debug version of the libraries (./configure --with-debugging=1) checks if the
942: index set is actually a permutation. The optimized version just believes you.
944: Note: The IS will be considered a permutation permanently, even if indices have been changes (for example, with
945: ISGeneralSetIndices()). It's a good idea to only set this property if the IS will not change in the future.
946: To clear this property, use ISClearInfoCache().
948: .seealso: ISPermutation(), ISSetInfo(), ISClearInfoCache().
949: @*/
950: PetscErrorCode ISSetPermutation(IS is)
951: {
953: if (PetscDefined(USE_DEBUG)) {
954: PetscMPIInt size;
956: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
957: if (size == 1) {
958: PetscInt i,n,*idx;
959: const PetscInt *iidx;
961: ISGetSize(is,&n);
962: PetscMalloc1(n,&idx);
963: ISGetIndices(is,&iidx);
964: PetscArraycpy(idx,iidx,n);
965: PetscIntSortSemiOrdered(n,idx);
966: for (i=0; i<n; i++) {
968: }
969: PetscFree(idx);
970: ISRestoreIndices(is,&iidx);
971: }
972: }
973: ISSetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_TRUE,PETSC_TRUE);
974: return 0;
975: }
977: /*@C
978: ISDestroy - Destroys an index set.
980: Collective on IS
982: Input Parameters:
983: . is - the index set
985: Level: beginner
987: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlocked()
988: @*/
989: PetscErrorCode ISDestroy(IS *is)
990: {
991: if (!*is) return 0;
993: if (--((PetscObject)(*is))->refct > 0) {*is = NULL; return 0;}
994: if ((*is)->complement) {
995: PetscInt refcnt;
996: PetscObjectGetReference((PetscObject)((*is)->complement), &refcnt);
998: ISDestroy(&(*is)->complement);
999: }
1000: if ((*is)->ops->destroy) {
1001: (*(*is)->ops->destroy)(*is);
1002: }
1003: PetscLayoutDestroy(&(*is)->map);
1004: /* Destroy local representations of offproc data. */
1005: PetscFree((*is)->total);
1006: PetscFree((*is)->nonlocal);
1007: PetscHeaderDestroy(is);
1008: return 0;
1009: }
1011: /*@
1012: ISInvertPermutation - Creates a new permutation that is the inverse of
1013: a given permutation.
1015: Collective on IS
1017: Input Parameters:
1018: + is - the index set
1019: - nlocal - number of indices on this processor in result (ignored for 1 proccessor) or
1020: use PETSC_DECIDE
1022: Output Parameter:
1023: . isout - the inverse permutation
1025: Level: intermediate
1027: Notes:
1028: For parallel index sets this does the complete parallel permutation, but the
1029: code is not efficient for huge index sets (10,000,000 indices).
1031: @*/
1032: PetscErrorCode ISInvertPermutation(IS is,PetscInt nlocal,IS *isout)
1033: {
1034: PetscBool isperm, isidentity, issame;
1038: ISGetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_TRUE,&isperm);
1040: ISGetInfo(is,IS_IDENTITY,IS_GLOBAL,PETSC_TRUE,&isidentity);
1041: issame = PETSC_FALSE;
1042: if (isidentity) {
1043: PetscInt n;
1044: PetscBool isallsame;
1046: ISGetLocalSize(is, &n);
1047: issame = (PetscBool) (n == nlocal);
1048: MPI_Allreduce(&issame, &isallsame, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is));
1049: issame = isallsame;
1050: }
1051: if (issame) {
1052: ISDuplicate(is,isout);
1053: } else {
1054: (*is->ops->invertpermutation)(is,nlocal,isout);
1055: ISSetPermutation(*isout);
1056: }
1057: return 0;
1058: }
1060: /*@
1061: ISGetSize - Returns the global length of an index set.
1063: Not Collective
1065: Input Parameter:
1066: . is - the index set
1068: Output Parameter:
1069: . size - the global size
1071: Level: beginner
1073: @*/
1074: PetscErrorCode ISGetSize(IS is,PetscInt *size)
1075: {
1078: *size = is->map->N;
1079: return 0;
1080: }
1082: /*@
1083: ISGetLocalSize - Returns the local (processor) length of an index set.
1085: Not Collective
1087: Input Parameter:
1088: . is - the index set
1090: Output Parameter:
1091: . size - the local size
1093: Level: beginner
1095: @*/
1096: PetscErrorCode ISGetLocalSize(IS is,PetscInt *size)
1097: {
1100: *size = is->map->n;
1101: return 0;
1102: }
1104: /*@
1105: ISGetLayout - get PetscLayout describing index set layout
1107: Not Collective
1109: Input Parameter:
1110: . is - the index set
1112: Output Parameter:
1113: . map - the layout
1115: Level: developer
1117: .seealso: ISSetLayout(), ISGetSize(), ISGetLocalSize()
1118: @*/
1119: PetscErrorCode ISGetLayout(IS is,PetscLayout *map)
1120: {
1123: *map = is->map;
1124: return 0;
1125: }
1127: /*@
1128: ISSetLayout - set PetscLayout describing index set layout
1130: Collective
1132: Input Arguments:
1133: + is - the index set
1134: - map - the layout
1136: Level: developer
1138: Notes:
1139: Users should typically use higher level functions such as ISCreateGeneral().
1141: This function can be useful in some special cases of constructing a new IS, e.g. after ISCreate() and before ISLoad().
1142: Otherwise, it is only valid to replace the layout with a layout known to be equivalent.
1144: .seealso: ISCreate(), ISGetLayout(), ISGetSize(), ISGetLocalSize()
1145: @*/
1146: PetscErrorCode ISSetLayout(IS is,PetscLayout map)
1147: {
1150: PetscLayoutReference(map,&is->map);
1151: return 0;
1152: }
1154: /*@C
1155: ISGetIndices - Returns a pointer to the indices. The user should call
1156: ISRestoreIndices() after having looked at the indices. The user should
1157: NOT change the indices.
1159: Not Collective
1161: Input Parameter:
1162: . is - the index set
1164: Output Parameter:
1165: . ptr - the location to put the pointer to the indices
1167: Fortran Note:
1168: This routine has two different interfaces from Fortran; the first is not recommend, it does not require Fortran 90
1169: $ IS is
1170: $ integer is_array(1)
1171: $ PetscOffset i_is
1172: $ int ierr
1173: $ call ISGetIndices(is,is_array,i_is,ierr)
1174: $
1175: $ Access first local entry in list
1176: $ value = is_array(i_is + 1)
1177: $
1178: $ ...... other code
1179: $ call ISRestoreIndices(is,is_array,i_is,ierr)
1180: The second Fortran interface is recommended.
1181: $ use petscisdef
1182: $ PetscInt, pointer :: array(:)
1183: $ PetscErrorCode ierr
1184: $ IS i
1185: $ 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
1192: .seealso: ISRestoreIndices(), ISGetIndicesF90()
1193: @*/
1194: PetscErrorCode ISGetIndices(IS is,const PetscInt *ptr[])
1195: {
1198: (*is->ops->getindices)(is,ptr);
1199: return 0;
1200: }
1202: /*@C
1203: ISGetMinMax - Gets the minimum and maximum values in an IS
1205: Not Collective
1207: Input Parameter:
1208: . is - the index set
1210: Output Parameters:
1211: + min - the minimum value
1212: - max - the maximum value
1214: Level: intermediate
1216: Notes:
1217: Empty index sets return min=PETSC_MAX_INT and max=PETSC_MIN_INT.
1218: In parallel, it returns the min and max of the local portion of the IS
1220: .seealso: ISGetIndices(), ISRestoreIndices(), ISGetIndicesF90()
1221: @*/
1222: PetscErrorCode ISGetMinMax(IS is,PetscInt *min,PetscInt *max)
1223: {
1225: if (min) *min = is->min;
1226: if (max) *max = is->max;
1227: return 0;
1228: }
1230: /*@
1231: ISLocate - determine the location of an index within the local component of an index set
1233: Not Collective
1235: Input Parameters:
1236: + is - the index set
1237: - key - the search key
1239: Output Parameter:
1240: . location - if >= 0, a location within the index set that is equal to the key, otherwise the key is not in the index set
1242: Level: intermediate
1243: @*/
1244: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
1245: {
1246: if (is->ops->locate) {
1247: (*is->ops->locate)(is,key,location);
1248: } else {
1249: PetscInt numIdx;
1250: PetscBool sorted;
1251: const PetscInt *idx;
1253: ISGetLocalSize(is,&numIdx);
1254: ISGetIndices(is,&idx);
1255: ISSorted(is,&sorted);
1256: if (sorted) {
1257: PetscFindInt(key,numIdx,idx,location);
1258: } else {
1259: PetscInt i;
1261: *location = -1;
1262: for (i = 0; i < numIdx; i++) {
1263: if (idx[i] == key) {
1264: *location = i;
1265: break;
1266: }
1267: }
1268: }
1269: ISRestoreIndices(is,&idx);
1270: }
1271: return 0;
1272: }
1274: /*@C
1275: ISRestoreIndices - Restores an index set to a usable state after a call
1276: to ISGetIndices().
1278: Not Collective
1280: Input Parameters:
1281: + is - the index set
1282: - ptr - the pointer obtained by ISGetIndices()
1284: Fortran Note:
1285: This routine is used differently from Fortran
1286: $ IS is
1287: $ integer is_array(1)
1288: $ PetscOffset i_is
1289: $ int ierr
1290: $ call ISGetIndices(is,is_array,i_is,ierr)
1291: $
1292: $ Access first local entry in list
1293: $ value = is_array(i_is + 1)
1294: $
1295: $ ...... other code
1296: $ call ISRestoreIndices(is,is_array,i_is,ierr)
1298: See the Fortran chapter of the users manual and
1299: petsc/src/vec/is/tests for details.
1301: Level: intermediate
1303: Note:
1304: This routine zeros out ptr. This is to prevent accidental us of the array after it has been restored.
1306: .seealso: ISGetIndices(), ISRestoreIndicesF90()
1307: @*/
1308: PetscErrorCode ISRestoreIndices(IS is,const PetscInt *ptr[])
1309: {
1312: if (is->ops->restoreindices) {
1313: (*is->ops->restoreindices)(is,ptr);
1314: }
1315: return 0;
1316: }
1318: static PetscErrorCode ISGatherTotal_Private(IS is)
1319: {
1320: PetscInt i,n,N;
1321: const PetscInt *lindices;
1322: MPI_Comm comm;
1323: PetscMPIInt rank,size,*sizes = NULL,*offsets = NULL,nn;
1327: PetscObjectGetComm((PetscObject)is,&comm);
1328: MPI_Comm_size(comm,&size);
1329: MPI_Comm_rank(comm,&rank);
1330: ISGetLocalSize(is,&n);
1331: PetscMalloc2(size,&sizes,size,&offsets);
1333: PetscMPIIntCast(n,&nn);
1334: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
1335: offsets[0] = 0;
1336: for (i=1; i<size; ++i) offsets[i] = offsets[i-1] + sizes[i-1];
1337: N = offsets[size-1] + sizes[size-1];
1339: PetscMalloc1(N,&(is->total));
1340: ISGetIndices(is,&lindices);
1341: MPI_Allgatherv((void*)lindices,nn,MPIU_INT,is->total,sizes,offsets,MPIU_INT,comm);
1342: ISRestoreIndices(is,&lindices);
1343: is->local_offset = offsets[rank];
1344: PetscFree2(sizes,offsets);
1345: return 0;
1346: }
1348: /*@C
1349: ISGetTotalIndices - Retrieve an array containing all indices across the communicator.
1351: Collective on IS
1353: Input Parameter:
1354: . is - the index set
1356: Output Parameter:
1357: . indices - total indices with rank 0 indices first, and so on; total array size is
1358: the same as returned with ISGetSize().
1360: Level: intermediate
1362: Notes:
1363: this is potentially nonscalable, but depends on the size of the total index set
1364: and the size of the communicator. This may be feasible for index sets defined on
1365: subcommunicators, such that the set size does not grow with PETSC_WORLD_COMM.
1366: Note also that there is no way to tell where the local part of the indices starts
1367: (use ISGetIndices() and ISGetNonlocalIndices() to retrieve just the local and just
1368: the nonlocal part (complement), respectively).
1370: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices(), ISGetSize()
1371: @*/
1372: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
1373: {
1374: PetscMPIInt size;
1378: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
1379: if (size == 1) {
1380: (*is->ops->getindices)(is,indices);
1381: } else {
1382: if (!is->total) {
1383: ISGatherTotal_Private(is);
1384: }
1385: *indices = is->total;
1386: }
1387: return 0;
1388: }
1390: /*@C
1391: ISRestoreTotalIndices - Restore the index array obtained with ISGetTotalIndices().
1393: Not Collective.
1395: Input Parameters:
1396: + is - the index set
1397: - indices - index array; must be the array obtained with ISGetTotalIndices()
1399: Level: intermediate
1401: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices()
1402: @*/
1403: PetscErrorCode ISRestoreTotalIndices(IS is, const PetscInt *indices[])
1404: {
1405: PetscMPIInt size;
1409: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
1410: if (size == 1) {
1411: (*is->ops->restoreindices)(is,indices);
1412: } else {
1414: }
1415: return 0;
1416: }
1418: /*@C
1419: ISGetNonlocalIndices - Retrieve an array of indices from remote processors
1420: in this communicator.
1422: Collective on IS
1424: Input Parameter:
1425: . is - the index set
1427: Output Parameter:
1428: . indices - indices with rank 0 indices first, and so on, omitting
1429: the current rank. Total number of indices is the difference
1430: total and local, obtained with ISGetSize() and ISGetLocalSize(),
1431: respectively.
1433: Level: intermediate
1435: Notes:
1436: restore the indices using ISRestoreNonlocalIndices().
1437: The same scalability considerations as those for ISGetTotalIndices
1438: apply here.
1440: .seealso: ISGetTotalIndices(), ISRestoreNonlocalIndices(), ISGetSize(), ISGetLocalSize().
1441: @*/
1442: PetscErrorCode ISGetNonlocalIndices(IS is, const PetscInt *indices[])
1443: {
1444: PetscMPIInt size;
1445: PetscInt n, N;
1449: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
1450: if (size == 1) *indices = NULL;
1451: else {
1452: if (!is->total) {
1453: ISGatherTotal_Private(is);
1454: }
1455: ISGetLocalSize(is,&n);
1456: ISGetSize(is,&N);
1457: PetscMalloc1(N-n, &(is->nonlocal));
1458: PetscArraycpy(is->nonlocal, is->total, is->local_offset);
1459: PetscArraycpy(is->nonlocal+is->local_offset, is->total+is->local_offset+n,N - is->local_offset - n);
1460: *indices = is->nonlocal;
1461: }
1462: return 0;
1463: }
1465: /*@C
1466: ISRestoreNonlocalIndices - Restore the index array obtained with ISGetNonlocalIndices().
1468: Not Collective.
1470: Input Parameters:
1471: + is - the index set
1472: - indices - index array; must be the array obtained with ISGetNonlocalIndices()
1474: Level: intermediate
1476: .seealso: ISGetTotalIndices(), ISGetNonlocalIndices(), ISRestoreTotalIndices()
1477: @*/
1478: PetscErrorCode ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
1479: {
1483: return 0;
1484: }
1486: /*@
1487: ISGetNonlocalIS - Gather all nonlocal indices for this IS and present
1488: them as another sequential index set.
1490: Collective on IS
1492: Input Parameter:
1493: . is - the index set
1495: Output Parameter:
1496: . complement - sequential IS with indices identical to the result of
1497: ISGetNonlocalIndices()
1499: Level: intermediate
1501: Notes:
1502: complement represents the result of ISGetNonlocalIndices as an IS.
1503: Therefore scalability issues similar to ISGetNonlocalIndices apply.
1504: The resulting IS must be restored using ISRestoreNonlocalIS().
1506: .seealso: ISGetNonlocalIndices(), ISRestoreNonlocalIndices(), ISAllGather(), ISGetSize()
1507: @*/
1508: PetscErrorCode ISGetNonlocalIS(IS is, IS *complement)
1509: {
1512: /* Check if the complement exists already. */
1513: if (is->complement) {
1514: *complement = is->complement;
1515: PetscObjectReference((PetscObject)(is->complement));
1516: } else {
1517: PetscInt N, n;
1518: const PetscInt *idx;
1519: ISGetSize(is, &N);
1520: ISGetLocalSize(is,&n);
1521: ISGetNonlocalIndices(is, &idx);
1522: ISCreateGeneral(PETSC_COMM_SELF, N-n,idx, PETSC_USE_POINTER, &(is->complement));
1523: PetscObjectReference((PetscObject)is->complement);
1524: *complement = is->complement;
1525: }
1526: return 0;
1527: }
1529: /*@
1530: ISRestoreNonlocalIS - Restore the IS obtained with ISGetNonlocalIS().
1532: Not collective.
1534: Input Parameters:
1535: + is - the index set
1536: - complement - index set of is's nonlocal indices
1538: Level: intermediate
1540: .seealso: ISGetNonlocalIS(), ISGetNonlocalIndices(), ISRestoreNonlocalIndices()
1541: @*/
1542: PetscErrorCode ISRestoreNonlocalIS(IS is, IS *complement)
1543: {
1544: PetscInt refcnt;
1549: PetscObjectGetReference((PetscObject)(is->complement), &refcnt);
1551: PetscObjectDereference((PetscObject)(is->complement));
1552: return 0;
1553: }
1555: /*@C
1556: ISViewFromOptions - View from Options
1558: Collective on IS
1560: Input Parameters:
1561: + A - the index set
1562: . obj - Optional object
1563: - name - command line option
1565: Level: intermediate
1566: .seealso: IS, ISView, PetscObjectViewFromOptions(), ISCreate()
1567: @*/
1568: PetscErrorCode ISViewFromOptions(IS A,PetscObject obj,const char name[])
1569: {
1571: PetscObjectViewFromOptions((PetscObject)A,obj,name);
1572: return 0;
1573: }
1575: /*@C
1576: ISView - Displays an index set.
1578: Collective on IS
1580: Input Parameters:
1581: + is - the index set
1582: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
1584: Level: intermediate
1586: .seealso: PetscViewerASCIIOpen()
1587: @*/
1588: PetscErrorCode ISView(IS is,PetscViewer viewer)
1589: {
1591: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is),&viewer);
1595: PetscObjectPrintClassNamePrefixType((PetscObject)is,viewer);
1596: PetscLogEventBegin(IS_View,is,viewer,0,0);
1597: (*is->ops->view)(is,viewer);
1598: PetscLogEventEnd(IS_View,is,viewer,0,0);
1599: return 0;
1600: }
1602: /*@
1603: ISLoad - Loads a vector that has been stored in binary or HDF5 format with ISView().
1605: Collective on PetscViewer
1607: Input Parameters:
1608: + is - the newly loaded vector, this needs to have been created with ISCreate() or some related function before a call to ISLoad().
1609: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or HDF5 file viewer, obtained from PetscViewerHDF5Open()
1611: Level: intermediate
1613: Notes:
1614: IF using HDF5, you must assign the IS the same name as was used in the IS
1615: that was stored in the file using PetscObjectSetName(). Otherwise you will
1616: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"
1618: .seealso: PetscViewerBinaryOpen(), ISView(), MatLoad(), VecLoad()
1619: @*/
1620: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
1621: {
1622: PetscBool isbinary, ishdf5;
1627: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
1628: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
1630: if (!((PetscObject)is)->type_name) ISSetType(is, ISGENERAL);
1631: PetscLogEventBegin(IS_Load,is,viewer,0,0);
1632: (*is->ops->load)(is, viewer);
1633: PetscLogEventEnd(IS_Load,is,viewer,0,0);
1634: return 0;
1635: }
1637: /*@
1638: ISSort - Sorts the indices of an index set.
1640: Collective on IS
1642: Input Parameters:
1643: . is - the index set
1645: Level: intermediate
1647: .seealso: ISSortRemoveDups(), ISSorted()
1648: @*/
1649: PetscErrorCode ISSort(IS is)
1650: {
1652: (*is->ops->sort)(is);
1653: ISSetInfo(is,IS_SORTED,IS_LOCAL,is->info_permanent[IS_LOCAL][IS_SORTED],PETSC_TRUE);
1654: return 0;
1655: }
1657: /*@
1658: ISSortRemoveDups - Sorts the indices of an index set, removing duplicates.
1660: Collective on IS
1662: Input Parameters:
1663: . is - the index set
1665: Level: intermediate
1667: .seealso: ISSort(), ISSorted()
1668: @*/
1669: PetscErrorCode ISSortRemoveDups(IS is)
1670: {
1672: ISClearInfoCache(is,PETSC_FALSE);
1673: (*is->ops->sortremovedups)(is);
1674: ISSetInfo(is,IS_SORTED,IS_LOCAL,is->info_permanent[IS_LOCAL][IS_SORTED],PETSC_TRUE);
1675: ISSetInfo(is,IS_UNIQUE,IS_LOCAL,is->info_permanent[IS_LOCAL][IS_UNIQUE],PETSC_TRUE);
1676: return 0;
1677: }
1679: /*@
1680: ISToGeneral - Converts an IS object of any type to ISGENERAL type
1682: Collective on IS
1684: Input Parameters:
1685: . is - the index set
1687: Level: intermediate
1689: .seealso: ISSorted()
1690: @*/
1691: PetscErrorCode ISToGeneral(IS is)
1692: {
1694: if (is->ops->togeneral) {
1695: (*is->ops->togeneral)(is);
1696: } else SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Not written for this type %s",((PetscObject)is)->type_name);
1697: return 0;
1698: }
1700: /*@
1701: ISSorted - Checks the indices to determine whether they have been sorted.
1703: Collective on IS
1705: Input Parameter:
1706: . is - the index set
1708: Output Parameter:
1709: . flg - output flag, either PETSC_TRUE if the index set is sorted,
1710: or PETSC_FALSE otherwise.
1712: Notes:
1713: For parallel IS objects this only indicates if the local part of the IS
1714: is sorted. So some processors may return PETSC_TRUE while others may
1715: return PETSC_FALSE.
1717: Level: intermediate
1719: .seealso: ISSort(), ISSortRemoveDups()
1720: @*/
1721: PetscErrorCode ISSorted(IS is,PetscBool *flg)
1722: {
1725: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,flg);
1726: return 0;
1727: }
1729: /*@
1730: ISDuplicate - Creates a duplicate copy of an index set.
1732: Collective on IS
1734: Input Parameter:
1735: . is - the index set
1737: Output Parameter:
1738: . isnew - the copy of the index set
1740: Level: beginner
1742: .seealso: ISCreateGeneral(), ISCopy()
1743: @*/
1744: PetscErrorCode ISDuplicate(IS is,IS *newIS)
1745: {
1748: (*is->ops->duplicate)(is,newIS);
1749: ISCopyInfo(is,*newIS);
1750: return 0;
1751: }
1753: /*@
1754: ISCopy - Copies an index set.
1756: Collective on IS
1758: Input Parameter:
1759: . is - the index set
1761: Output Parameter:
1762: . isy - the copy of the index set
1764: Level: beginner
1766: .seealso: ISDuplicate()
1767: @*/
1768: PetscErrorCode ISCopy(IS is,IS isy)
1769: {
1773: if (is == isy) return 0;
1774: ISCopyInfo(is,isy);
1775: isy->max = is->max;
1776: isy->min = is->min;
1777: (*is->ops->copy)(is,isy);
1778: return 0;
1779: }
1781: /*@
1782: ISOnComm - Split a parallel IS on subcomms (usually self) or concatenate index sets on subcomms into a parallel index set
1784: Collective on IS
1786: Input Parameters:
1787: + is - index set
1788: . comm - communicator for new index set
1789: - mode - copy semantics, PETSC_USE_POINTER for no-copy if possible, otherwise PETSC_COPY_VALUES
1791: Output Parameter:
1792: . newis - new IS on comm
1794: Level: advanced
1796: Notes:
1797: It is usually desirable to create a parallel IS and look at the local part when necessary.
1799: This function is useful if serial ISs must be created independently, or to view many
1800: logically independent serial ISs.
1802: The input IS must have the same type on every process.
1803: @*/
1804: PetscErrorCode ISOnComm(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
1805: {
1806: PetscMPIInt match;
1810: MPI_Comm_compare(PetscObjectComm((PetscObject)is),comm,&match);
1811: if (mode != PETSC_COPY_VALUES && (match == MPI_IDENT || match == MPI_CONGRUENT)) {
1812: PetscObjectReference((PetscObject)is);
1813: *newis = is;
1814: } else {
1815: (*is->ops->oncomm)(is,comm,mode,newis);
1816: }
1817: return 0;
1818: }
1820: /*@
1821: ISSetBlockSize - informs an index set that it has a given block size
1823: Logicall Collective on IS
1825: Input Parameters:
1826: + is - index set
1827: - bs - block size
1829: Level: intermediate
1831: Notes:
1832: This is much like the block size for Vecs. It indicates that one can think of the indices as
1833: being in a collection of equal size blocks. For ISBlock() these collections of blocks are all contiquous
1834: within a block but this is not the case for other IS.
1835: ISBlockGetIndices() only works for ISBlock IS, not others.
1837: .seealso: ISGetBlockSize(), ISCreateBlock(), ISBlockGetIndices(),
1838: @*/
1839: PetscErrorCode ISSetBlockSize(IS is,PetscInt bs)
1840: {
1844: if (PetscDefined(USE_DEBUG)) {
1845: const PetscInt *indices;
1846: PetscInt length,i,j;
1847: ISGetIndices(is,&indices);
1848: if (indices) {
1849: ISGetLocalSize(is,&length);
1851: for (i=0;i<length/bs;i+=bs) {
1852: for (j=0;j<bs-1;j++) {
1854: }
1855: }
1856: }
1857: ISRestoreIndices(is,&indices);
1858: }
1859: (*is->ops->setblocksize)(is,bs);
1860: return 0;
1861: }
1863: /*@
1864: ISGetBlockSize - Returns the number of elements in a block.
1866: Not Collective
1868: Input Parameter:
1869: . is - the index set
1871: Output Parameter:
1872: . size - the number of elements in a block
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: ISBlockGetSize(), ISGetSize(), ISCreateBlock(), ISSetBlockSize()
1883: @*/
1884: PetscErrorCode ISGetBlockSize(IS is,PetscInt *size)
1885: {
1886: PetscLayoutGetBlockSize(is->map, size);
1887: return 0;
1888: }
1890: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[])
1891: {
1892: PetscInt len,i;
1893: const PetscInt *ptr;
1895: ISGetLocalSize(is,&len);
1896: ISGetIndices(is,&ptr);
1897: for (i=0; i<len; i++) idx[i] = ptr[i];
1898: ISRestoreIndices(is,&ptr);
1899: return 0;
1900: }
1902: /*MC
1903: ISGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1904: The users should call ISRestoreIndicesF90() after having looked at the
1905: indices. The user should NOT change the indices.
1907: Synopsis:
1908: ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1910: Not collective
1912: Input Parameter:
1913: . x - index set
1915: Output Parameters:
1916: + xx_v - the Fortran90 pointer to the array
1917: - ierr - error code
1919: Example of Usage:
1920: .vb
1921: PetscInt, pointer xx_v(:)
1922: ....
1923: call ISGetIndicesF90(x,xx_v,ierr)
1924: a = xx_v(3)
1925: call ISRestoreIndicesF90(x,xx_v,ierr)
1926: .ve
1928: Level: intermediate
1930: .seealso: ISRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices()
1932: M*/
1934: /*MC
1935: ISRestoreIndicesF90 - Restores an index set to a usable state after
1936: a call to ISGetIndicesF90().
1938: Synopsis:
1939: ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1941: Not collective
1943: Input Parameters:
1944: + x - index set
1945: - xx_v - the Fortran90 pointer to the array
1947: Output Parameter:
1948: . ierr - error code
1950: Example of Usage:
1951: .vb
1952: PetscInt, pointer xx_v(:)
1953: ....
1954: call ISGetIndicesF90(x,xx_v,ierr)
1955: a = xx_v(3)
1956: call ISRestoreIndicesF90(x,xx_v,ierr)
1957: .ve
1959: Level: intermediate
1961: .seealso: ISGetIndicesF90(), ISGetIndices(), ISRestoreIndices()
1963: M*/
1965: /*MC
1966: ISBlockGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1967: The users should call ISBlockRestoreIndicesF90() after having looked at the
1968: indices. The user should NOT change the indices.
1970: Synopsis:
1971: ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1973: Not collective
1975: Input Parameter:
1976: . x - index set
1978: Output Parameters:
1979: + xx_v - the Fortran90 pointer to the array
1980: - ierr - error code
1981: Example of Usage:
1982: .vb
1983: PetscInt, pointer xx_v(:)
1984: ....
1985: call ISBlockGetIndicesF90(x,xx_v,ierr)
1986: a = xx_v(3)
1987: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
1988: .ve
1990: Level: intermediate
1992: .seealso: ISBlockRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices(),
1993: ISRestoreIndices()
1995: M*/
1997: /*MC
1998: ISBlockRestoreIndicesF90 - Restores an index set to a usable state after
1999: a call to ISBlockGetIndicesF90().
2001: Synopsis:
2002: ISBlockRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
2004: Not Collective
2006: Input Parameters:
2007: + x - index set
2008: - xx_v - the Fortran90 pointer to the array
2010: Output Parameter:
2011: . ierr - error code
2013: Example of Usage:
2014: .vb
2015: PetscInt, pointer xx_v(:)
2016: ....
2017: call ISBlockGetIndicesF90(x,xx_v,ierr)
2018: a = xx_v(3)
2019: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
2020: .ve
2022: Notes:
2023: Not yet supported for all F90 compilers
2025: Level: intermediate
2027: .seealso: ISBlockGetIndicesF90(), ISGetIndices(), ISRestoreIndices(), ISRestoreIndicesF90()
2029: M*/