Actual source code: general.c
1: #define PETSCVEC_DLL
2: /*
3: Provides the functions for index sets (IS) defined by a list of integers.
4: */
5: #include src/vec/is/impls/general/general.h
6: #include petscvec.h
10: PetscErrorCode ISDuplicate_General(IS is,IS *newIS)
11: {
13: IS_General *sub = (IS_General *)is->data;
16: ISCreateGeneral(is->comm,sub->n,sub->idx,newIS);
17: return(0);
18: }
22: PetscErrorCode ISDestroy_General(IS is)
23: {
24: IS_General *is_general = (IS_General*)is->data;
28: if (is_general->allocated) {
29: PetscFree(is_general->idx);
30: }
31: PetscFree(is_general);
32: PetscHeaderDestroy(is);
33: return(0);
34: }
38: PetscErrorCode ISIdentity_General(IS is,PetscTruth *ident)
39: {
40: IS_General *is_general = (IS_General*)is->data;
41: PetscInt i,n = is_general->n,*idx = is_general->idx;
44: is->isidentity = PETSC_TRUE;
45: *ident = PETSC_TRUE;
46: for (i=0; i<n; i++) {
47: if (idx[i] != i) {
48: is->isidentity = PETSC_FALSE;
49: *ident = PETSC_FALSE;
50: break;
51: }
52: }
53: return(0);
54: }
58: PetscErrorCode ISGetIndices_General(IS in,PetscInt **idx)
59: {
60: IS_General *sub = (IS_General*)in->data;
63: *idx = sub->idx;
64: return(0);
65: }
69: PetscErrorCode ISRestoreIndices_General(IS in,PetscInt **idx)
70: {
71: IS_General *sub = (IS_General*)in->data;
74: if (*idx != sub->idx) {
75: SETERRQ(PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
76: }
77: return(0);
78: }
82: PetscErrorCode ISGetSize_General(IS is,PetscInt *size)
83: {
84: IS_General *sub = (IS_General *)is->data;
87: *size = sub->N;
88: return(0);
89: }
93: PetscErrorCode ISGetLocalSize_General(IS is,PetscInt *size)
94: {
95: IS_General *sub = (IS_General *)is->data;
98: *size = sub->n;
99: return(0);
100: }
104: PetscErrorCode ISInvertPermutation_General(IS is,PetscInt nlocal,IS *isout)
105: {
106: IS_General *sub = (IS_General *)is->data;
107: PetscInt i,*ii,n = sub->n,*idx = sub->idx,nstart;
108: PetscMPIInt size;
109: IS istmp,nistmp;
113: MPI_Comm_size(is->comm,&size);
114: if (size == 1) {
115: PetscMalloc(n*sizeof(PetscInt),&ii);
116: for (i=0; i<n; i++) {
117: ii[idx[i]] = i;
118: }
119: ISCreateGeneral(PETSC_COMM_SELF,n,ii,isout);
120: ISSetPermutation(*isout);
121: PetscFree(ii);
122: } else {
123: /* crude, nonscalable get entire IS on each processor */
124: if (nlocal == PETSC_DECIDE) SETERRQ(PETSC_ERR_SUP,"Do not yet support nlocal of PETSC_DECIDE");
125: ISAllGather(is,&istmp);
126: ISSetPermutation(istmp);
127: ISInvertPermutation(istmp,PETSC_DECIDE,&nistmp);
128: ISDestroy(istmp);
129: /* get the part we need */
130: MPI_Scan(&nlocal,&nstart,1,MPIU_INT,MPI_SUM,is->comm);
131: #if defined(PETSC_USE_DEBUG)
132: {
133: PetscMPIInt rank;
134: MPI_Comm_rank(is->comm,&rank);
135: if (rank == size-1) {
136: if (nstart != sub->N) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Sum of nlocal lengths %d != total IS length %d",nstart,sub->N);
137: }
138: }
139: #endif
140: nstart -= nlocal;
141: ISGetIndices(nistmp,&idx);
142: ISCreateGeneral(is->comm,nlocal,idx+nstart,isout);
143: ISRestoreIndices(nistmp,&idx);
144: ISDestroy(nistmp);
145: }
146: return(0);
147: }
151: PetscErrorCode ISView_General(IS is,PetscViewer viewer)
152: {
153: IS_General *sub = (IS_General *)is->data;
155: PetscInt i,n = sub->n,*idx = sub->idx;
156: PetscTruth iascii;
159: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
160: if (iascii) {
161: MPI_Comm comm;
162: PetscMPIInt rank,size;
164: PetscObjectGetComm((PetscObject)viewer,&comm);
165: MPI_Comm_rank(comm,&rank);
166: MPI_Comm_size(comm,&size);
168: if (size > 1) {
169: if (is->isperm) {
170: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutation\n",rank);
171: }
172: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in set %D\n",rank,n);
173: for (i=0; i<n; i++) {
174: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,idx[i]);
175: }
176: } else {
177: if (is->isperm) {
178: PetscViewerASCIISynchronizedPrintf(viewer,"Index set is permutation\n");
179: }
180: PetscViewerASCIISynchronizedPrintf(viewer,"Number of indices in set %D\n",n);
181: for (i=0; i<n; i++) {
182: PetscViewerASCIISynchronizedPrintf(viewer,"%D %D\n",i,idx[i]);
183: }
184: }
185: PetscViewerFlush(viewer);
186: } else {
187: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for this object",((PetscObject)viewer)->type_name);
188: }
189: return(0);
190: }
194: PetscErrorCode ISSort_General(IS is)
195: {
196: IS_General *sub = (IS_General *)is->data;
200: if (sub->sorted) return(0);
201: PetscSortInt(sub->n,sub->idx);
202: sub->sorted = PETSC_TRUE;
203: return(0);
204: }
208: PetscErrorCode ISSorted_General(IS is,PetscTruth *flg)
209: {
210: IS_General *sub = (IS_General *)is->data;
213: *flg = sub->sorted;
214: return(0);
215: }
217: static struct _ISOps myops = { ISGetSize_General,
218: ISGetLocalSize_General,
219: ISGetIndices_General,
220: ISRestoreIndices_General,
221: ISInvertPermutation_General,
222: ISSort_General,
223: ISSorted_General,
224: ISDuplicate_General,
225: ISDestroy_General,
226: ISView_General,
227: ISIdentity_General };
231: PetscErrorCode ISCreateGeneral_Private(MPI_Comm comm,IS *is)
232: {
234: IS Nindex = *is;
235: IS_General *sub = (IS_General*)Nindex->data;
236: PetscInt n = sub->n,i,min,max;
237: const PetscInt *idx = sub->idx;
238: PetscTruth sorted = PETSC_TRUE;
239: PetscTruth flg;
243: if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
245: *is = PETSC_NULL;
246: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
247: VecInitializePackage(PETSC_NULL);
248: #endif
250: MPI_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,comm);
251: for (i=1; i<n; i++) {
252: if (idx[i] < idx[i-1]) {sorted = PETSC_FALSE; break;}
253: }
254: if (n) {min = max = idx[0];} else {min = max = 0;}
255: for (i=1; i<n; i++) {
256: if (idx[i] < min) min = idx[i];
257: if (idx[i] > max) max = idx[i];
258: }
259: sub->sorted = sorted;
260: Nindex->min = min;
261: Nindex->max = max;
262: PetscMemcpy(Nindex->ops,&myops,sizeof(myops));
263: Nindex->isperm = PETSC_FALSE;
264: Nindex->isidentity = PETSC_FALSE;
265: PetscOptionsHasName(PETSC_NULL,"-is_view",&flg);
266: if (flg) {
267: PetscViewer viewer;
268: PetscViewerASCIIGetStdout(Nindex->comm,&viewer);
269: ISView(Nindex,viewer);
270: }
271: *is = Nindex;
272: return(0);
273: }
275: /*@
276: ISCreateGeneral - Creates a data structure for an index set
277: containing a list of integers.
279: Collective on MPI_Comm
281: Input Parameters:
282: + comm - the MPI communicator
283: . n - the length of the index set
284: - idx - the list of integers
286: Output Parameter:
287: . is - the new index set
289: Notes:
290: The index array is copied to internally allocated storage. After the call,
291: the user can free the index array.
293: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
294: conceptually the same as MPI_Group operations. The IS are then
295: distributed sets of indices and thus certain operations on them are
296: collective.
298: Level: beginner
300: Concepts: index sets^creating
301: Concepts: IS^creating
303: .seealso: ISCreateGeneralWithArray(), ISCreateStride(), ISCreateBlock(), ISAllGather()
304: @*/
305: PetscErrorCode ISCreateGeneral(MPI_Comm comm,PetscInt n,const PetscInt idx[],IS *is)
306: {
308: IS Nindex;
309: IS_General *sub;
313: if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
315: *is = PETSC_NULL;
316: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
317: VecInitializePackage(PETSC_NULL);
318: #endif
320: PetscHeaderCreate(Nindex,_p_IS,struct _ISOps,IS_COOKIE,IS_GENERAL,"IS",comm,ISDestroy,ISView);
321: PetscNew(IS_General,&sub);
322: PetscLogObjectMemory(Nindex,sizeof(IS_General)+n*sizeof(PetscInt)+sizeof(struct _p_IS));
323: PetscMalloc(n*sizeof(PetscInt),&sub->idx);
324: PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));
325: sub->n = n;
326: sub->allocated = PETSC_TRUE;
327: Nindex->data = (void*)sub;
329: *is = Nindex;
330: ISCreateGeneral_Private(comm,is);
332: return(0);
333: }
335: /*@C
336: ISCreateGeneralWithArray - Creates a data structure for an index set
337: containing a list of integers.
339: Collective on MPI_Comm
341: Input Parameters:
342: + comm - the MPI communicator
343: . n - the length of the index set
344: - idx - the list of integers
346: Output Parameter:
347: . is - the new index set
349: Notes:
350: Unlike with ISCreateGeneral, the indices are not copied to internally
351: allocated storage. The user array is not freed by ISDestroy.
353: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
354: conceptually the same as MPI_Group operations. The IS are then
355: distributed sets of indices and thus certain operations on them are collective.
357: Level: beginner
359: Concepts: index sets^creating
360: Concepts: IS^creating
362: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
363: @*/
364: PetscErrorCode ISCreateGeneralWithArray(MPI_Comm comm,PetscInt n,PetscInt idx[],IS *is)
365: {
367: IS Nindex;
368: IS_General *sub;
372: if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
374: *is = PETSC_NULL;
375: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
376: VecInitializePackage(PETSC_NULL);
377: #endif
379: PetscHeaderCreate(Nindex,_p_IS,struct _ISOps,IS_COOKIE,IS_GENERAL,"IS",comm,ISDestroy,ISView);
380: PetscNew(IS_General,&sub);
381: PetscLogObjectMemory(Nindex,sizeof(IS_General)+n*sizeof(PetscInt)+sizeof(struct _p_IS));
382: sub->idx = idx;
383: sub->n = n;
384: sub->allocated = PETSC_FALSE;
385: Nindex->data = (void*)sub;
387: *is = Nindex;
388: ISCreateGeneral_Private(comm,is);
390: return(0);
391: }