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: ISGetIndices(is, &idx);
647: for (i = 1; i < n; i++) if (idx[i] != idx[i-1] + 1) break;
648: if (i < n) intervalLocal = PETSC_FALSE;
649: ISRestoreIndices(is, &idx);
650: }
652: if (type == IS_LOCAL || size == 1) {
653: *flg = intervalLocal;
654: } else {
655: MPI_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
656: if (*flg) {
657: PetscInt n, min = PETSC_MAX_INT, max = PETSC_MIN_INT;
658: PetscInt maxprev;
660: ISGetLocalSize(is, &n);
661: if (n) ISGetMinMax(is, &min, &max);
662: maxprev = PETSC_MIN_INT;
663: MPI_Exscan(&max, &maxprev, 1, MPIU_INT, MPI_MAX, comm);
664: if (rank && n && (maxprev != min - 1)) intervalLocal = PETSC_FALSE;
665: MPI_Allreduce(&intervalLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
666: }
667: }
668: }
669: return 0;
670: }
672: static PetscErrorCode ISGetInfo_Identity(IS is, ISInfoType type, PetscBool *flg)
673: {
674: MPI_Comm comm;
675: PetscMPIInt size, rank;
677: comm = PetscObjectComm((PetscObject)is);
678: MPI_Comm_size(comm, &size);
679: MPI_Comm_size(comm, &rank);
680: if (type == IS_GLOBAL && is->ops->intervalglobal) {
681: PetscBool isinterval;
683: (*is->ops->intervalglobal)(is,&isinterval);
684: *flg = PETSC_FALSE;
685: if (isinterval) {
686: PetscInt min;
688: ISGetMinMax(is, &min, NULL);
689: MPI_Bcast(&min, 1, MPIU_INT, 0, comm);
690: if (min == 0) *flg = PETSC_TRUE;
691: }
692: } else if (type == IS_LOCAL && is->ops->intervallocal) {
693: PetscBool isinterval;
695: (*is->ops->intervallocal)(is,&isinterval);
696: *flg = PETSC_FALSE;
697: if (isinterval) {
698: PetscInt min;
700: ISGetMinMax(is, &min, NULL);
701: if (min == 0) *flg = PETSC_TRUE;
702: }
703: } else {
704: PetscBool identLocal;
705: PetscInt n, i, rStart;
706: const PetscInt *idx;
708: ISGetLocalSize(is, &n);
709: ISGetIndices(is, &idx);
710: PetscLayoutGetRange(is->map, &rStart, NULL);
711: identLocal = PETSC_TRUE;
712: for (i = 0; i < n; i++) {
713: if (idx[i] != rStart + i) break;
714: }
715: if (i < n) identLocal = PETSC_FALSE;
716: if (type == IS_LOCAL || size == 1) {
717: *flg = identLocal;
718: } else {
719: MPI_Allreduce(&identLocal, flg, 1, MPIU_BOOL, MPI_LAND, comm);
720: }
721: ISRestoreIndices(is, &idx);
722: }
723: return 0;
724: }
726: /*@
727: ISGetInfo - Determine whether an index set satisfies a given property
729: 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())
731: Input Parameters:
732: + is - the index set
733: . info - describing a property of the index set, one of those listed in the documentation of ISSetInfo()
734: . 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
735: - type - whether the property is local (IS_LOCAL) or global (IS_GLOBAL)
737: Output Parameter:
738: . flg - wheter the property is true (PETSC_TRUE) or false (PETSC_FALSE)
740: 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().
742: Level: advanced
744: .seealso: ISInfo, ISInfoType, ISSetInfo(), ISClearInfoCache()
746: @*/
747: PetscErrorCode ISGetInfo(IS is, ISInfo info, ISInfoType type, PetscBool compute, PetscBool *flg)
748: {
749: MPI_Comm comm, errcomm;
750: PetscMPIInt rank, size;
751: PetscInt itype;
752: PetscBool hasprop;
753: PetscBool infer;
757: comm = PetscObjectComm((PetscObject)is);
758: if (type == IS_GLOBAL) {
760: errcomm = comm;
761: } else {
762: errcomm = PETSC_COMM_SELF;
763: }
765: MPI_Comm_size(comm, &size);
766: MPI_Comm_rank(comm, &rank);
769: if (size == 1) type = IS_LOCAL;
770: itype = (type == IS_LOCAL) ? 0 : 1;
771: hasprop = PETSC_FALSE;
772: infer = PETSC_FALSE;
773: if (is->info_permanent[itype][(int)info]) {
774: hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
775: infer = PETSC_TRUE;
776: } else if ((itype == IS_LOCAL) && (is->info[IS_LOCAL][info] != IS_INFO_UNKNOWN)) {
777: /* we can cache local properties as long as we clear them when the IS changes */
778: /* NOTE: we only cache local values because there is no ISAssemblyBegin()/ISAssemblyEnd(),
779: so we have no way of knowing when a cached value has been invalidated by changes on a different process */
780: hasprop = (is->info[itype][(int)info] == IS_INFO_TRUE) ? PETSC_TRUE : PETSC_FALSE;
781: infer = PETSC_TRUE;
782: } else if (compute) {
783: switch (info) {
784: case IS_SORTED:
785: ISGetInfo_Sorted(is, type, &hasprop);
786: break;
787: case IS_UNIQUE:
788: ISGetInfo_Unique(is, type, &hasprop);
789: break;
790: case IS_PERMUTATION:
791: ISGetInfo_Permutation(is, type, &hasprop);
792: break;
793: case IS_INTERVAL:
794: ISGetInfo_Interval(is, type, &hasprop);
795: break;
796: case IS_IDENTITY:
797: ISGetInfo_Identity(is, type, &hasprop);
798: break;
799: default:
800: SETERRQ(errcomm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown IS property");
801: }
802: infer = PETSC_TRUE;
803: }
804: /* call ISSetInfo_Internal to keep all of the implications straight */
805: if (infer) ISSetInfo_Internal(is, info, type, IS_INFO_UNKNOWN, hasprop);
806: *flg = hasprop;
807: return 0;
808: }
810: static PetscErrorCode ISCopyInfo(IS source, IS dest)
811: {
812: PetscArraycpy(&dest->info[0], &source->info[0], 2);
813: PetscArraycpy(&dest->info_permanent[0], &source->info_permanent[0], 2);
814: return 0;
815: }
817: /*@
818: ISIdentity - Determines whether index set is the identity mapping.
820: Collective on IS
822: Input Parameters:
823: . is - the index set
825: Output Parameters:
826: . ident - PETSC_TRUE if an identity, else PETSC_FALSE
828: Level: intermediate
830: Note: If ISSetIdentity() (or ISSetInfo() for a permanent property) has been called,
831: ISIdentity() will return its answer without communication between processes, but
832: otherwise the output ident will be computed from ISGetInfo(),
833: which may require synchronization on the communicator of IS. To avoid this computation,
834: call ISGetInfo() directly with the compute flag set to PETSC_FALSE, and ident will be assumed false.
836: .seealso: ISSetIdentity(), ISGetInfo()
837: @*/
838: PetscErrorCode ISIdentity(IS is,PetscBool *ident)
839: {
842: ISGetInfo(is,IS_IDENTITY,IS_GLOBAL,PETSC_TRUE,ident);
843: return 0;
844: }
846: /*@
847: ISSetIdentity - Informs the index set that it is an identity.
849: Logically Collective on IS
851: Input Parameter:
852: . is - the index set
854: Level: intermediate
856: Note: The IS will be considered the identity permanently, even if indices have been changes (for example, with
857: ISGeneralSetIndices()). It's a good idea to only set this property if the IS will not change in the future.
858: To clear this property, use ISClearInfoCache().
860: .seealso: ISIdentity(), ISSetInfo(), ISClearInfoCache()
861: @*/
862: PetscErrorCode ISSetIdentity(IS is)
863: {
865: ISSetInfo(is,IS_IDENTITY,IS_GLOBAL,PETSC_TRUE,PETSC_TRUE);
866: return 0;
867: }
869: /*@
870: ISContiguousLocal - Locates an index set with contiguous range within a global range, if possible
872: Not Collective
874: Input Parameters:
875: + is - the index set
876: . gstart - global start
877: - gend - global end
879: Output Parameters:
880: + start - start of contiguous block, as an offset from gstart
881: - contig - PETSC_TRUE if the index set refers to contiguous entries on this process, else PETSC_FALSE
883: Level: developer
885: .seealso: ISGetLocalSize(), VecGetOwnershipRange()
886: @*/
887: PetscErrorCode ISContiguousLocal(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
888: {
892: *start = -1;
893: *contig = PETSC_FALSE;
894: if (is->ops->contiguous) {
895: (*is->ops->contiguous)(is,gstart,gend,start,contig);
896: }
897: return 0;
898: }
900: /*@
901: ISPermutation - PETSC_TRUE or PETSC_FALSE depending on whether the
902: index set has been declared to be a permutation.
904: Logically Collective on IS
906: Input Parameter:
907: . is - the index set
909: Output Parameter:
910: . perm - PETSC_TRUE if a permutation, else PETSC_FALSE
912: Level: intermediate
914: Note: If it is not alread known that the IS is a permutation (if ISSetPermutation()
915: or ISSetInfo() has not been called), this routine will not attempt to compute
916: whether the index set is a permutation and will assume perm is PETSC_FALSE.
917: To compute the value when it is not already known, use ISGetInfo() with
918: the compute flag set to PETSC_TRUE.
920: .seealso: ISSetPermutation(), ISGetInfo()
921: @*/
922: PetscErrorCode ISPermutation(IS is,PetscBool *perm)
923: {
926: ISGetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_FALSE,perm);
927: return 0;
928: }
930: /*@
931: ISSetPermutation - Informs the index set that it is a permutation.
933: Logically Collective on IS
935: Input Parameter:
936: . is - the index set
938: Level: intermediate
940: The debug version of the libraries (./configure --with-debugging=1) checks if the
941: index set is actually a permutation. The optimized version just believes you.
943: Note: The IS will be considered a permutation permanently, even if indices have been changes (for example, with
944: ISGeneralSetIndices()). It's a good idea to only set this property if the IS will not change in the future.
945: To clear this property, use ISClearInfoCache().
947: .seealso: ISPermutation(), ISSetInfo(), ISClearInfoCache().
948: @*/
949: PetscErrorCode ISSetPermutation(IS is)
950: {
952: if (PetscDefined(USE_DEBUG)) {
953: PetscMPIInt size;
955: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
956: if (size == 1) {
957: PetscInt i,n,*idx;
958: const PetscInt *iidx;
960: ISGetSize(is,&n);
961: PetscMalloc1(n,&idx);
962: ISGetIndices(is,&iidx);
963: PetscArraycpy(idx,iidx,n);
964: PetscIntSortSemiOrdered(n,idx);
965: for (i=0; i<n; i++) {
967: }
968: PetscFree(idx);
969: ISRestoreIndices(is,&iidx);
970: }
971: }
972: ISSetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_TRUE,PETSC_TRUE);
973: return 0;
974: }
976: /*@C
977: ISDestroy - Destroys an index set.
979: Collective on IS
981: Input Parameters:
982: . is - the index set
984: Level: beginner
986: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlocked()
987: @*/
988: PetscErrorCode ISDestroy(IS *is)
989: {
990: if (!*is) return 0;
992: if (--((PetscObject)(*is))->refct > 0) {*is = NULL; return 0;}
993: if ((*is)->complement) {
994: PetscInt refcnt;
995: PetscObjectGetReference((PetscObject)((*is)->complement), &refcnt);
997: ISDestroy(&(*is)->complement);
998: }
999: if ((*is)->ops->destroy) {
1000: (*(*is)->ops->destroy)(*is);
1001: }
1002: PetscLayoutDestroy(&(*is)->map);
1003: /* Destroy local representations of offproc data. */
1004: PetscFree((*is)->total);
1005: PetscFree((*is)->nonlocal);
1006: PetscHeaderDestroy(is);
1007: return 0;
1008: }
1010: /*@
1011: ISInvertPermutation - Creates a new permutation that is the inverse of
1012: a given permutation.
1014: Collective on IS
1016: Input Parameters:
1017: + is - the index set
1018: - nlocal - number of indices on this processor in result (ignored for 1 proccessor) or
1019: use PETSC_DECIDE
1021: Output Parameter:
1022: . isout - the inverse permutation
1024: Level: intermediate
1026: Notes:
1027: For parallel index sets this does the complete parallel permutation, but the
1028: code is not efficient for huge index sets (10,000,000 indices).
1030: @*/
1031: PetscErrorCode ISInvertPermutation(IS is,PetscInt nlocal,IS *isout)
1032: {
1033: PetscBool isperm, isidentity, issame;
1037: ISGetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_TRUE,&isperm);
1039: ISGetInfo(is,IS_IDENTITY,IS_GLOBAL,PETSC_TRUE,&isidentity);
1040: issame = PETSC_FALSE;
1041: if (isidentity) {
1042: PetscInt n;
1043: PetscBool isallsame;
1045: ISGetLocalSize(is, &n);
1046: issame = (PetscBool) (n == nlocal);
1047: MPI_Allreduce(&issame, &isallsame, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)is));
1048: issame = isallsame;
1049: }
1050: if (issame) {
1051: ISDuplicate(is,isout);
1052: } else {
1053: (*is->ops->invertpermutation)(is,nlocal,isout);
1054: ISSetPermutation(*isout);
1055: }
1056: return 0;
1057: }
1059: /*@
1060: ISGetSize - Returns the global length of an index set.
1062: Not Collective
1064: Input Parameter:
1065: . is - the index set
1067: Output Parameter:
1068: . size - the global size
1070: Level: beginner
1072: @*/
1073: PetscErrorCode ISGetSize(IS is,PetscInt *size)
1074: {
1077: *size = is->map->N;
1078: return 0;
1079: }
1081: /*@
1082: ISGetLocalSize - Returns the local (processor) length of an index set.
1084: Not Collective
1086: Input Parameter:
1087: . is - the index set
1089: Output Parameter:
1090: . size - the local size
1092: Level: beginner
1094: @*/
1095: PetscErrorCode ISGetLocalSize(IS is,PetscInt *size)
1096: {
1099: *size = is->map->n;
1100: return 0;
1101: }
1103: /*@
1104: ISGetLayout - get PetscLayout describing index set layout
1106: Not Collective
1108: Input Parameter:
1109: . is - the index set
1111: Output Parameter:
1112: . map - the layout
1114: Level: developer
1116: .seealso: ISSetLayout(), ISGetSize(), ISGetLocalSize()
1117: @*/
1118: PetscErrorCode ISGetLayout(IS is,PetscLayout *map)
1119: {
1122: *map = is->map;
1123: return 0;
1124: }
1126: /*@
1127: ISSetLayout - set PetscLayout describing index set layout
1129: Collective
1131: Input Arguments:
1132: + is - the index set
1133: - map - the layout
1135: Level: developer
1137: Notes:
1138: Users should typically use higher level functions such as ISCreateGeneral().
1140: This function can be useful in some special cases of constructing a new IS, e.g. after ISCreate() and before ISLoad().
1141: Otherwise, it is only valid to replace the layout with a layout known to be equivalent.
1143: .seealso: ISCreate(), ISGetLayout(), ISGetSize(), ISGetLocalSize()
1144: @*/
1145: PetscErrorCode ISSetLayout(IS is,PetscLayout map)
1146: {
1149: PetscLayoutReference(map,&is->map);
1150: return 0;
1151: }
1153: /*@C
1154: ISGetIndices - Returns a pointer to the indices. The user should call
1155: ISRestoreIndices() after having looked at the indices. The user should
1156: NOT change the indices.
1158: Not Collective
1160: Input Parameter:
1161: . is - the index set
1163: Output Parameter:
1164: . ptr - the location to put the pointer to the indices
1166: Fortran Note:
1167: This routine has two different interfaces from Fortran; the first is not recommend, it does not require Fortran 90
1168: $ IS is
1169: $ integer is_array(1)
1170: $ PetscOffset i_is
1171: $ int ierr
1172: $ call ISGetIndices(is,is_array,i_is,ierr)
1173: $
1174: $ Access first local entry in list
1175: $ value = is_array(i_is + 1)
1176: $
1177: $ ...... other code
1178: $ call ISRestoreIndices(is,is_array,i_is,ierr)
1179: The second Fortran interface is recommended.
1180: $ use petscisdef
1181: $ PetscInt, pointer :: array(:)
1182: $ PetscErrorCode ierr
1183: $ IS i
1184: $ call ISGetIndicesF90(i,array,ierr)
1186: See the Fortran chapter of the users manual and
1187: petsc/src/is/[tutorials,tests] for details.
1189: Level: intermediate
1191: .seealso: ISRestoreIndices(), ISGetIndicesF90()
1192: @*/
1193: PetscErrorCode ISGetIndices(IS is,const PetscInt *ptr[])
1194: {
1197: (*is->ops->getindices)(is,ptr);
1198: return 0;
1199: }
1201: /*@C
1202: ISGetMinMax - Gets the minimum and maximum values in an IS
1204: Not Collective
1206: Input Parameter:
1207: . is - the index set
1209: Output Parameters:
1210: + min - the minimum value
1211: - max - the maximum value
1213: Level: intermediate
1215: Notes:
1216: Empty index sets return min=PETSC_MAX_INT and max=PETSC_MIN_INT.
1217: In parallel, it returns the min and max of the local portion of the IS
1219: .seealso: ISGetIndices(), ISRestoreIndices(), ISGetIndicesF90()
1220: @*/
1221: PetscErrorCode ISGetMinMax(IS is,PetscInt *min,PetscInt *max)
1222: {
1224: if (min) *min = is->min;
1225: if (max) *max = is->max;
1226: return 0;
1227: }
1229: /*@
1230: ISLocate - determine the location of an index within the local component of an index set
1232: Not Collective
1234: Input Parameters:
1235: + is - the index set
1236: - key - the search key
1238: Output Parameter:
1239: . location - if >= 0, a location within the index set that is equal to the key, otherwise the key is not in the index set
1241: Level: intermediate
1242: @*/
1243: PetscErrorCode ISLocate(IS is, PetscInt key, PetscInt *location)
1244: {
1245: if (is->ops->locate) {
1246: (*is->ops->locate)(is,key,location);
1247: } else {
1248: PetscInt numIdx;
1249: PetscBool sorted;
1250: const PetscInt *idx;
1252: ISGetLocalSize(is,&numIdx);
1253: ISGetIndices(is,&idx);
1254: ISSorted(is,&sorted);
1255: if (sorted) {
1256: PetscFindInt(key,numIdx,idx,location);
1257: } else {
1258: PetscInt i;
1260: *location = -1;
1261: for (i = 0; i < numIdx; i++) {
1262: if (idx[i] == key) {
1263: *location = i;
1264: break;
1265: }
1266: }
1267: }
1268: ISRestoreIndices(is,&idx);
1269: }
1270: return 0;
1271: }
1273: /*@C
1274: ISRestoreIndices - Restores an index set to a usable state after a call
1275: to ISGetIndices().
1277: Not Collective
1279: Input Parameters:
1280: + is - the index set
1281: - ptr - the pointer obtained by ISGetIndices()
1283: Fortran Note:
1284: This routine is used differently from Fortran
1285: $ IS is
1286: $ integer is_array(1)
1287: $ PetscOffset i_is
1288: $ int ierr
1289: $ call ISGetIndices(is,is_array,i_is,ierr)
1290: $
1291: $ Access first local entry in list
1292: $ value = is_array(i_is + 1)
1293: $
1294: $ ...... other code
1295: $ call ISRestoreIndices(is,is_array,i_is,ierr)
1297: See the Fortran chapter of the users manual and
1298: petsc/src/vec/is/tests for details.
1300: Level: intermediate
1302: Note:
1303: This routine zeros out ptr. This is to prevent accidental us of the array after it has been restored.
1305: .seealso: ISGetIndices(), ISRestoreIndicesF90()
1306: @*/
1307: PetscErrorCode ISRestoreIndices(IS is,const PetscInt *ptr[])
1308: {
1311: if (is->ops->restoreindices) {
1312: (*is->ops->restoreindices)(is,ptr);
1313: }
1314: return 0;
1315: }
1317: static PetscErrorCode ISGatherTotal_Private(IS is)
1318: {
1319: PetscInt i,n,N;
1320: const PetscInt *lindices;
1321: MPI_Comm comm;
1322: PetscMPIInt rank,size,*sizes = NULL,*offsets = NULL,nn;
1326: PetscObjectGetComm((PetscObject)is,&comm);
1327: MPI_Comm_size(comm,&size);
1328: MPI_Comm_rank(comm,&rank);
1329: ISGetLocalSize(is,&n);
1330: PetscMalloc2(size,&sizes,size,&offsets);
1332: PetscMPIIntCast(n,&nn);
1333: MPI_Allgather(&nn,1,MPI_INT,sizes,1,MPI_INT,comm);
1334: offsets[0] = 0;
1335: for (i=1; i<size; ++i) offsets[i] = offsets[i-1] + sizes[i-1];
1336: N = offsets[size-1] + sizes[size-1];
1338: PetscMalloc1(N,&(is->total));
1339: ISGetIndices(is,&lindices);
1340: MPI_Allgatherv((void*)lindices,nn,MPIU_INT,is->total,sizes,offsets,MPIU_INT,comm);
1341: ISRestoreIndices(is,&lindices);
1342: is->local_offset = offsets[rank];
1343: PetscFree2(sizes,offsets);
1344: return 0;
1345: }
1347: /*@C
1348: ISGetTotalIndices - Retrieve an array containing all indices across the communicator.
1350: Collective on IS
1352: Input Parameter:
1353: . is - the index set
1355: Output Parameter:
1356: . indices - total indices with rank 0 indices first, and so on; total array size is
1357: the same as returned with ISGetSize().
1359: Level: intermediate
1361: Notes:
1362: this is potentially nonscalable, but depends on the size of the total index set
1363: and the size of the communicator. This may be feasible for index sets defined on
1364: subcommunicators, such that the set size does not grow with PETSC_WORLD_COMM.
1365: Note also that there is no way to tell where the local part of the indices starts
1366: (use ISGetIndices() and ISGetNonlocalIndices() to retrieve just the local and just
1367: the nonlocal part (complement), respectively).
1369: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices(), ISGetSize()
1370: @*/
1371: PetscErrorCode ISGetTotalIndices(IS is, const PetscInt *indices[])
1372: {
1373: PetscMPIInt size;
1377: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
1378: if (size == 1) {
1379: (*is->ops->getindices)(is,indices);
1380: } else {
1381: if (!is->total) {
1382: ISGatherTotal_Private(is);
1383: }
1384: *indices = is->total;
1385: }
1386: return 0;
1387: }
1389: /*@C
1390: ISRestoreTotalIndices - Restore the index array obtained with ISGetTotalIndices().
1392: Not Collective.
1394: Input Parameters:
1395: + is - the index set
1396: - indices - index array; must be the array obtained with ISGetTotalIndices()
1398: Level: intermediate
1400: .seealso: ISRestoreTotalIndices(), ISGetNonlocalIndices()
1401: @*/
1402: PetscErrorCode ISRestoreTotalIndices(IS is, const PetscInt *indices[])
1403: {
1404: PetscMPIInt size;
1408: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
1409: if (size == 1) {
1410: (*is->ops->restoreindices)(is,indices);
1411: } else {
1413: }
1414: return 0;
1415: }
1417: /*@C
1418: ISGetNonlocalIndices - Retrieve an array of indices from remote processors
1419: in this communicator.
1421: Collective on IS
1423: Input Parameter:
1424: . is - the index set
1426: Output Parameter:
1427: . indices - indices with rank 0 indices first, and so on, omitting
1428: the current rank. Total number of indices is the difference
1429: total and local, obtained with ISGetSize() and ISGetLocalSize(),
1430: respectively.
1432: Level: intermediate
1434: Notes:
1435: restore the indices using ISRestoreNonlocalIndices().
1436: The same scalability considerations as those for ISGetTotalIndices
1437: apply here.
1439: .seealso: ISGetTotalIndices(), ISRestoreNonlocalIndices(), ISGetSize(), ISGetLocalSize().
1440: @*/
1441: PetscErrorCode ISGetNonlocalIndices(IS is, const PetscInt *indices[])
1442: {
1443: PetscMPIInt size;
1444: PetscInt n, N;
1448: MPI_Comm_size(PetscObjectComm((PetscObject)is), &size);
1449: if (size == 1) *indices = NULL;
1450: else {
1451: if (!is->total) {
1452: ISGatherTotal_Private(is);
1453: }
1454: ISGetLocalSize(is,&n);
1455: ISGetSize(is,&N);
1456: PetscMalloc1(N-n, &(is->nonlocal));
1457: PetscArraycpy(is->nonlocal, is->total, is->local_offset);
1458: PetscArraycpy(is->nonlocal+is->local_offset, is->total+is->local_offset+n,N - is->local_offset - n);
1459: *indices = is->nonlocal;
1460: }
1461: return 0;
1462: }
1464: /*@C
1465: ISRestoreNonlocalIndices - Restore the index array obtained with ISGetNonlocalIndices().
1467: Not Collective.
1469: Input Parameters:
1470: + is - the index set
1471: - indices - index array; must be the array obtained with ISGetNonlocalIndices()
1473: Level: intermediate
1475: .seealso: ISGetTotalIndices(), ISGetNonlocalIndices(), ISRestoreTotalIndices()
1476: @*/
1477: PetscErrorCode ISRestoreNonlocalIndices(IS is, const PetscInt *indices[])
1478: {
1482: return 0;
1483: }
1485: /*@
1486: ISGetNonlocalIS - Gather all nonlocal indices for this IS and present
1487: them as another sequential index set.
1489: Collective on IS
1491: Input Parameter:
1492: . is - the index set
1494: Output Parameter:
1495: . complement - sequential IS with indices identical to the result of
1496: ISGetNonlocalIndices()
1498: Level: intermediate
1500: Notes:
1501: complement represents the result of ISGetNonlocalIndices as an IS.
1502: Therefore scalability issues similar to ISGetNonlocalIndices apply.
1503: The resulting IS must be restored using ISRestoreNonlocalIS().
1505: .seealso: ISGetNonlocalIndices(), ISRestoreNonlocalIndices(), ISAllGather(), ISGetSize()
1506: @*/
1507: PetscErrorCode ISGetNonlocalIS(IS is, IS *complement)
1508: {
1511: /* Check if the complement exists already. */
1512: if (is->complement) {
1513: *complement = is->complement;
1514: PetscObjectReference((PetscObject)(is->complement));
1515: } else {
1516: PetscInt N, n;
1517: const PetscInt *idx;
1518: ISGetSize(is, &N);
1519: ISGetLocalSize(is,&n);
1520: ISGetNonlocalIndices(is, &idx);
1521: ISCreateGeneral(PETSC_COMM_SELF, N-n,idx, PETSC_USE_POINTER, &(is->complement));
1522: PetscObjectReference((PetscObject)is->complement);
1523: *complement = is->complement;
1524: }
1525: return 0;
1526: }
1528: /*@
1529: ISRestoreNonlocalIS - Restore the IS obtained with ISGetNonlocalIS().
1531: Not collective.
1533: Input Parameters:
1534: + is - the index set
1535: - complement - index set of is's nonlocal indices
1537: Level: intermediate
1539: .seealso: ISGetNonlocalIS(), ISGetNonlocalIndices(), ISRestoreNonlocalIndices()
1540: @*/
1541: PetscErrorCode ISRestoreNonlocalIS(IS is, IS *complement)
1542: {
1543: PetscInt refcnt;
1548: PetscObjectGetReference((PetscObject)(is->complement), &refcnt);
1550: PetscObjectDereference((PetscObject)(is->complement));
1551: return 0;
1552: }
1554: /*@C
1555: ISViewFromOptions - View from Options
1557: Collective on IS
1559: Input Parameters:
1560: + A - the index set
1561: . obj - Optional object
1562: - name - command line option
1564: Level: intermediate
1565: .seealso: IS, ISView, PetscObjectViewFromOptions(), ISCreate()
1566: @*/
1567: PetscErrorCode ISViewFromOptions(IS A,PetscObject obj,const char name[])
1568: {
1570: PetscObjectViewFromOptions((PetscObject)A,obj,name);
1571: return 0;
1572: }
1574: /*@C
1575: ISView - Displays an index set.
1577: Collective on IS
1579: Input Parameters:
1580: + is - the index set
1581: - viewer - viewer used to display the set, for example PETSC_VIEWER_STDOUT_SELF.
1583: Level: intermediate
1585: .seealso: PetscViewerASCIIOpen()
1586: @*/
1587: PetscErrorCode ISView(IS is,PetscViewer viewer)
1588: {
1590: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is),&viewer);
1594: PetscObjectPrintClassNamePrefixType((PetscObject)is,viewer);
1595: PetscLogEventBegin(IS_View,is,viewer,0,0);
1596: (*is->ops->view)(is,viewer);
1597: PetscLogEventEnd(IS_View,is,viewer,0,0);
1598: return 0;
1599: }
1601: /*@
1602: ISLoad - Loads a vector that has been stored in binary or HDF5 format with ISView().
1604: Collective on PetscViewer
1606: Input Parameters:
1607: + is - the newly loaded vector, this needs to have been created with ISCreate() or some related function before a call to ISLoad().
1608: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or HDF5 file viewer, obtained from PetscViewerHDF5Open()
1610: Level: intermediate
1612: Notes:
1613: IF using HDF5, you must assign the IS the same name as was used in the IS
1614: that was stored in the file using PetscObjectSetName(). Otherwise you will
1615: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"
1617: .seealso: PetscViewerBinaryOpen(), ISView(), MatLoad(), VecLoad()
1618: @*/
1619: PetscErrorCode ISLoad(IS is, PetscViewer viewer)
1620: {
1621: PetscBool isbinary, ishdf5;
1626: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
1627: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
1629: if (!((PetscObject)is)->type_name) ISSetType(is, ISGENERAL);
1630: PetscLogEventBegin(IS_Load,is,viewer,0,0);
1631: (*is->ops->load)(is, viewer);
1632: PetscLogEventEnd(IS_Load,is,viewer,0,0);
1633: return 0;
1634: }
1636: /*@
1637: ISSort - Sorts the indices of an index set.
1639: Collective on IS
1641: Input Parameters:
1642: . is - the index set
1644: Level: intermediate
1646: .seealso: ISSortRemoveDups(), ISSorted()
1647: @*/
1648: PetscErrorCode ISSort(IS is)
1649: {
1651: (*is->ops->sort)(is);
1652: ISSetInfo(is,IS_SORTED,IS_LOCAL,is->info_permanent[IS_LOCAL][IS_SORTED],PETSC_TRUE);
1653: return 0;
1654: }
1656: /*@
1657: ISSortRemoveDups - Sorts the indices of an index set, removing duplicates.
1659: Collective on IS
1661: Input Parameters:
1662: . is - the index set
1664: Level: intermediate
1666: .seealso: ISSort(), ISSorted()
1667: @*/
1668: PetscErrorCode ISSortRemoveDups(IS is)
1669: {
1671: ISClearInfoCache(is,PETSC_FALSE);
1672: (*is->ops->sortremovedups)(is);
1673: ISSetInfo(is,IS_SORTED,IS_LOCAL,is->info_permanent[IS_LOCAL][IS_SORTED],PETSC_TRUE);
1674: ISSetInfo(is,IS_UNIQUE,IS_LOCAL,is->info_permanent[IS_LOCAL][IS_UNIQUE],PETSC_TRUE);
1675: return 0;
1676: }
1678: /*@
1679: ISToGeneral - Converts an IS object of any type to ISGENERAL type
1681: Collective on IS
1683: Input Parameters:
1684: . is - the index set
1686: Level: intermediate
1688: .seealso: ISSorted()
1689: @*/
1690: PetscErrorCode ISToGeneral(IS is)
1691: {
1693: if (is->ops->togeneral) {
1694: (*is->ops->togeneral)(is);
1695: } else SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_SUP,"Not written for this type %s",((PetscObject)is)->type_name);
1696: return 0;
1697: }
1699: /*@
1700: ISSorted - Checks the indices to determine whether they have been sorted.
1702: Collective on IS
1704: Input Parameter:
1705: . is - the index set
1707: Output Parameter:
1708: . flg - output flag, either PETSC_TRUE if the index set is sorted,
1709: or PETSC_FALSE otherwise.
1711: Notes:
1712: For parallel IS objects this only indicates if the local part of the IS
1713: is sorted. So some processors may return PETSC_TRUE while others may
1714: return PETSC_FALSE.
1716: Level: intermediate
1718: .seealso: ISSort(), ISSortRemoveDups()
1719: @*/
1720: PetscErrorCode ISSorted(IS is,PetscBool *flg)
1721: {
1724: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,flg);
1725: return 0;
1726: }
1728: /*@
1729: ISDuplicate - Creates a duplicate copy of an index set.
1731: Collective on IS
1733: Input Parameter:
1734: . is - the index set
1736: Output Parameter:
1737: . isnew - the copy of the index set
1739: Level: beginner
1741: .seealso: ISCreateGeneral(), ISCopy()
1742: @*/
1743: PetscErrorCode ISDuplicate(IS is,IS *newIS)
1744: {
1747: (*is->ops->duplicate)(is,newIS);
1748: ISCopyInfo(is,*newIS);
1749: return 0;
1750: }
1752: /*@
1753: ISCopy - Copies an index set.
1755: Collective on IS
1757: Input Parameter:
1758: . is - the index set
1760: Output Parameter:
1761: . isy - the copy of the index set
1763: Level: beginner
1765: .seealso: ISDuplicate()
1766: @*/
1767: PetscErrorCode ISCopy(IS is,IS isy)
1768: {
1772: if (is == isy) return 0;
1773: ISCopyInfo(is,isy);
1774: isy->max = is->max;
1775: isy->min = is->min;
1776: (*is->ops->copy)(is,isy);
1777: return 0;
1778: }
1780: /*@
1781: ISOnComm - Split a parallel IS on subcomms (usually self) or concatenate index sets on subcomms into a parallel index set
1783: Collective on IS
1785: Input Parameters:
1786: + is - index set
1787: . comm - communicator for new index set
1788: - mode - copy semantics, PETSC_USE_POINTER for no-copy if possible, otherwise PETSC_COPY_VALUES
1790: Output Parameter:
1791: . newis - new IS on comm
1793: Level: advanced
1795: Notes:
1796: It is usually desirable to create a parallel IS and look at the local part when necessary.
1798: This function is useful if serial ISs must be created independently, or to view many
1799: logically independent serial ISs.
1801: The input IS must have the same type on every process.
1802: @*/
1803: PetscErrorCode ISOnComm(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
1804: {
1805: PetscMPIInt match;
1809: MPI_Comm_compare(PetscObjectComm((PetscObject)is),comm,&match);
1810: if (mode != PETSC_COPY_VALUES && (match == MPI_IDENT || match == MPI_CONGRUENT)) {
1811: PetscObjectReference((PetscObject)is);
1812: *newis = is;
1813: } else {
1814: (*is->ops->oncomm)(is,comm,mode,newis);
1815: }
1816: return 0;
1817: }
1819: /*@
1820: ISSetBlockSize - informs an index set that it has a given block size
1822: Logicall Collective on IS
1824: Input Parameters:
1825: + is - index set
1826: - bs - block size
1828: Level: intermediate
1830: Notes:
1831: This is much like the block size for Vecs. It indicates that one can think of the indices as
1832: being in a collection of equal size blocks. For ISBlock() these collections of blocks are all contiquous
1833: within a block but this is not the case for other IS.
1834: ISBlockGetIndices() only works for ISBlock IS, not others.
1836: .seealso: ISGetBlockSize(), ISCreateBlock(), ISBlockGetIndices(),
1837: @*/
1838: PetscErrorCode ISSetBlockSize(IS is,PetscInt bs)
1839: {
1843: if (PetscDefined(USE_DEBUG)) {
1844: const PetscInt *indices;
1845: PetscInt length,i,j;
1846: ISGetIndices(is,&indices);
1847: if (indices) {
1848: ISGetLocalSize(is,&length);
1850: for (i=0;i<length/bs;i+=bs) {
1851: for (j=0;j<bs-1;j++) {
1853: }
1854: }
1855: }
1856: ISRestoreIndices(is,&indices);
1857: }
1858: (*is->ops->setblocksize)(is,bs);
1859: return 0;
1860: }
1862: /*@
1863: ISGetBlockSize - Returns the number of elements in a block.
1865: Not Collective
1867: Input Parameter:
1868: . is - the index set
1870: Output Parameter:
1871: . size - the number of elements in a block
1873: Level: intermediate
1875: Notes:
1876: This is much like the block size for Vecs. It indicates that one can think of the indices as
1877: being in a collection of equal size blocks. For ISBlock() these collections of blocks are all contiquous
1878: within a block but this is not the case for other IS.
1879: ISBlockGetIndices() only works for ISBlock IS, not others.
1881: .seealso: ISBlockGetSize(), ISGetSize(), ISCreateBlock(), ISSetBlockSize()
1882: @*/
1883: PetscErrorCode ISGetBlockSize(IS is,PetscInt *size)
1884: {
1885: PetscLayoutGetBlockSize(is->map, size);
1886: return 0;
1887: }
1889: PetscErrorCode ISGetIndicesCopy(IS is, PetscInt idx[])
1890: {
1891: PetscInt len,i;
1892: const PetscInt *ptr;
1894: ISGetLocalSize(is,&len);
1895: ISGetIndices(is,&ptr);
1896: for (i=0; i<len; i++) idx[i] = ptr[i];
1897: ISRestoreIndices(is,&ptr);
1898: return 0;
1899: }
1901: /*MC
1902: ISGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1903: The users should call ISRestoreIndicesF90() after having looked at the
1904: indices. The user should NOT change the indices.
1906: Synopsis:
1907: ISGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1909: Not collective
1911: Input Parameter:
1912: . x - index set
1914: Output Parameters:
1915: + xx_v - the Fortran90 pointer to the array
1916: - ierr - error code
1918: Example of Usage:
1919: .vb
1920: PetscInt, pointer xx_v(:)
1921: ....
1922: call ISGetIndicesF90(x,xx_v,ierr)
1923: a = xx_v(3)
1924: call ISRestoreIndicesF90(x,xx_v,ierr)
1925: .ve
1927: Level: intermediate
1929: .seealso: ISRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices()
1931: M*/
1933: /*MC
1934: ISRestoreIndicesF90 - Restores an index set to a usable state after
1935: a call to ISGetIndicesF90().
1937: Synopsis:
1938: ISRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1940: Not collective
1942: Input Parameters:
1943: + x - index set
1944: - xx_v - the Fortran90 pointer to the array
1946: Output Parameter:
1947: . ierr - error code
1949: Example of Usage:
1950: .vb
1951: PetscInt, pointer xx_v(:)
1952: ....
1953: call ISGetIndicesF90(x,xx_v,ierr)
1954: a = xx_v(3)
1955: call ISRestoreIndicesF90(x,xx_v,ierr)
1956: .ve
1958: Level: intermediate
1960: .seealso: ISGetIndicesF90(), ISGetIndices(), ISRestoreIndices()
1962: M*/
1964: /*MC
1965: ISBlockGetIndicesF90 - Accesses the elements of an index set from Fortran90.
1966: The users should call ISBlockRestoreIndicesF90() after having looked at the
1967: indices. The user should NOT change the indices.
1969: Synopsis:
1970: ISBlockGetIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
1972: Not collective
1974: Input Parameter:
1975: . x - index set
1977: Output Parameters:
1978: + xx_v - the Fortran90 pointer to the array
1979: - ierr - error code
1980: Example of Usage:
1981: .vb
1982: PetscInt, pointer xx_v(:)
1983: ....
1984: call ISBlockGetIndicesF90(x,xx_v,ierr)
1985: a = xx_v(3)
1986: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
1987: .ve
1989: Level: intermediate
1991: .seealso: ISBlockRestoreIndicesF90(), ISGetIndices(), ISRestoreIndices(),
1992: ISRestoreIndices()
1994: M*/
1996: /*MC
1997: ISBlockRestoreIndicesF90 - Restores an index set to a usable state after
1998: a call to ISBlockGetIndicesF90().
2000: Synopsis:
2001: ISBlockRestoreIndicesF90(IS x,{integer, pointer :: xx_v(:)},integer ierr)
2003: Not Collective
2005: Input Parameters:
2006: + x - index set
2007: - xx_v - the Fortran90 pointer to the array
2009: Output Parameter:
2010: . ierr - error code
2012: Example of Usage:
2013: .vb
2014: PetscInt, pointer xx_v(:)
2015: ....
2016: call ISBlockGetIndicesF90(x,xx_v,ierr)
2017: a = xx_v(3)
2018: call ISBlockRestoreIndicesF90(x,xx_v,ierr)
2019: .ve
2021: Notes:
2022: Not yet supported for all F90 compilers
2024: Level: intermediate
2026: .seealso: ISBlockGetIndicesF90(), ISGetIndices(), ISRestoreIndices(), ISRestoreIndicesF90()
2028: M*/