Actual source code: general.c
petsc-3.4.5 2014-06-29
2: /*
3: Provides the functions for index sets (IS) defined by a list of integers.
4: */
5: #include <../src/vec/is/is/impls/general/general.h> /*I "petscis.h" I*/
6: #include <petscvec.h>
7: #include <petscviewer.h>
11: PetscErrorCode ISDuplicate_General(IS is,IS *newIS)
12: {
14: IS_General *sub = (IS_General*)is->data;
17: ISCreateGeneral(PetscObjectComm((PetscObject)is),sub->n,sub->idx,PETSC_COPY_VALUES,newIS);
18: return(0);
19: }
23: PetscErrorCode ISDestroy_General(IS is)
24: {
25: IS_General *is_general = (IS_General*)is->data;
29: if (is_general->allocated) {
30: PetscFree(is_general->idx);
31: }
32: PetscObjectComposeFunction((PetscObject)is,"ISGeneralSetIndices_C",0);
33: PetscFree(is->data);
34: return(0);
35: }
39: PetscErrorCode ISIdentity_General(IS is,PetscBool *ident)
40: {
41: IS_General *is_general = (IS_General*)is->data;
42: PetscInt i,n = is_general->n,*idx = is_general->idx;
45: is->isidentity = PETSC_TRUE;
46: *ident = PETSC_TRUE;
47: for (i=0; i<n; i++) {
48: if (idx[i] != i) {
49: is->isidentity = PETSC_FALSE;
50: *ident = PETSC_FALSE;
51: break;
52: }
53: }
54: return(0);
55: }
59: static PetscErrorCode ISCopy_General(IS is,IS isy)
60: {
61: IS_General *is_general = (IS_General*)is->data,*isy_general = (IS_General*)isy->data;
65: if (is_general->n != isy_general->n || is_general->N != isy_general->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Index sets incompatible");
66: isy_general->sorted = is_general->sorted;
67: PetscMemcpy(isy_general->idx,is_general->idx,is_general->n*sizeof(PetscInt));
68: return(0);
69: }
73: PetscErrorCode ISOnComm_General(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
74: {
76: IS_General *sub = (IS_General*)is->data;
79: if (mode == PETSC_OWN_POINTER) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Cannot use PETSC_OWN_POINTER");
80: ISCreateGeneral(comm,sub->n,sub->idx,mode,newis);
81: return(0);
82: }
86: static PetscErrorCode ISSetBlockSize_General(IS is,PetscInt bs)
87: {
88: IS_General *sub = (IS_General*)is->data;
91: if (sub->N % bs) SETERRQ2(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_SIZ,"Block size %D does not divide global size %D",bs,sub->N);
92: if (sub->n % bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Block size %D does not divide local size %D",bs,sub->n);
93: #if defined(PETSC_USE_DEBUG)
94: {
95: PetscInt i,j;
96: for (i=0; i<sub->n; i+=bs) {
97: for (j=0; j<bs; j++) {
98: if (sub->idx[i+j] != sub->idx[i]+j) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Index set does not have block structure, cannot set block size to %D",bs);
99: }
100: }
101: }
102: #endif
103: is->bs = bs;
104: return(0);
105: }
109: static PetscErrorCode ISContiguousLocal_General(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
110: {
111: IS_General *sub = (IS_General*)is->data;
112: PetscInt i,p;
115: *start = 0;
116: *contig = PETSC_TRUE;
117: if (!sub->n) return(0);
118: p = sub->idx[0];
119: if (p < gstart) goto nomatch;
120: *start = p - gstart;
121: if (sub->n > gend-p) goto nomatch;
122: for (i=1; i<sub->n; i++,p++) {
123: if (sub->idx[i] != p+1) goto nomatch;
124: }
125: return(0);
126: nomatch:
127: *start = -1;
128: *contig = PETSC_FALSE;
129: return(0);
130: }
134: PetscErrorCode ISGetIndices_General(IS in,const PetscInt *idx[])
135: {
136: IS_General *sub = (IS_General*)in->data;
139: *idx = sub->idx;
140: return(0);
141: }
145: PetscErrorCode ISRestoreIndices_General(IS in,const PetscInt *idx[])
146: {
147: IS_General *sub = (IS_General*)in->data;
150: if (*idx != sub->idx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
151: return(0);
152: }
156: PetscErrorCode ISGetSize_General(IS is,PetscInt *size)
157: {
158: IS_General *sub = (IS_General*)is->data;
161: *size = sub->N;
162: return(0);
163: }
167: PetscErrorCode ISGetLocalSize_General(IS is,PetscInt *size)
168: {
169: IS_General *sub = (IS_General*)is->data;
172: *size = sub->n;
173: return(0);
174: }
178: PetscErrorCode ISInvertPermutation_General(IS is,PetscInt nlocal,IS *isout)
179: {
180: IS_General *sub = (IS_General*)is->data;
181: PetscInt i,*ii,n = sub->n,nstart;
182: const PetscInt *idx = sub->idx;
183: PetscMPIInt size;
184: IS istmp,nistmp;
188: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
189: if (size == 1) {
190: PetscMalloc(n*sizeof(PetscInt),&ii);
191: for (i=0; i<n; i++) ii[idx[i]] = i;
192: ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_OWN_POINTER,isout);
193: ISSetPermutation(*isout);
194: } else {
195: /* crude, nonscalable get entire IS on each processor */
196: if (nlocal == PETSC_DECIDE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Do not yet support nlocal of PETSC_DECIDE");
197: ISAllGather(is,&istmp);
198: ISSetPermutation(istmp);
199: ISInvertPermutation(istmp,PETSC_DECIDE,&nistmp);
200: ISDestroy(&istmp);
201: /* get the part we need */
202: MPI_Scan(&nlocal,&nstart,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)is));
203: #if defined(PETSC_USE_DEBUG)
204: {
205: PetscMPIInt rank;
206: MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
207: if (rank == size-1) {
208: if (nstart != sub->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Sum of nlocal lengths %d != total IS length %d",nstart,sub->N);
209: }
210: }
211: #endif
212: nstart -= nlocal;
213: ISGetIndices(nistmp,&idx);
214: ISCreateGeneral(PetscObjectComm((PetscObject)is),nlocal,idx+nstart,PETSC_COPY_VALUES,isout);
215: ISRestoreIndices(nistmp,&idx);
216: ISDestroy(&nistmp);
217: }
218: return(0);
219: }
223: PetscErrorCode ISView_General_Binary(IS is,PetscViewer viewer)
224: {
226: IS_General *isa = (IS_General*) is->data;
227: PetscMPIInt rank,size,mesgsize,tag = ((PetscObject)viewer)->tag, mesglen;
228: PetscInt len,j,tr[2];
229: int fdes;
230: MPI_Status status;
231: PetscInt message_count,flowcontrolcount,*values;
234: PetscViewerBinaryGetDescriptor(viewer,&fdes);
236: /* determine maximum message to arrive */
237: MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
238: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
240: tr[0] = IS_FILE_CLASSID;
241: tr[1] = isa->N;
242: PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_FALSE);
243: MPI_Reduce(&isa->n,&len,1,MPIU_INT,MPI_SUM,0,PetscObjectComm((PetscObject)is));
245: PetscViewerFlowControlStart(viewer,&message_count,&flowcontrolcount);
246: if (!rank) {
247: PetscBinaryWrite(fdes,isa->idx,isa->n,PETSC_INT,PETSC_FALSE);
249: PetscMalloc(len*sizeof(PetscInt),&values);
250: PetscMPIIntCast(len,&mesgsize);
251: /* receive and save messages */
252: for (j=1; j<size; j++) {
253: PetscViewerFlowControlStepMaster(viewer,j,&message_count,flowcontrolcount);
254: MPI_Recv(values,mesgsize,MPIU_INT,j,tag,PetscObjectComm((PetscObject)is),&status);
255: MPI_Get_count(&status,MPIU_INT,&mesglen);
256: PetscBinaryWrite(fdes,values,(PetscInt)mesglen,PETSC_INT,PETSC_TRUE);
257: }
258: PetscViewerFlowControlEndMaster(viewer,&message_count);
259: PetscFree(values);
260: } else {
261: PetscViewerFlowControlStepWorker(viewer,rank,&message_count);
262: PetscMPIIntCast(isa->n,&mesgsize);
263: MPI_Send(isa->idx,mesgsize,MPIU_INT,0,tag,PetscObjectComm((PetscObject)is));
264: PetscViewerFlowControlEndWorker(viewer,&message_count);
265: }
266: return(0);
267: }
271: PetscErrorCode ISView_General(IS is,PetscViewer viewer)
272: {
273: IS_General *sub = (IS_General*)is->data;
275: PetscInt i,n = sub->n,*idx = sub->idx;
276: PetscBool iascii,isbinary;
279: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
280: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
281: if (iascii) {
282: MPI_Comm comm;
283: PetscMPIInt rank,size;
285: PetscObjectGetComm((PetscObject)viewer,&comm);
286: MPI_Comm_rank(comm,&rank);
287: MPI_Comm_size(comm,&size);
289: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
290: if (size > 1) {
291: if (is->isperm) {
292: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutation\n",rank);
293: }
294: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in set %D\n",rank,n);
295: for (i=0; i<n; i++) {
296: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,idx[i]);
297: }
298: } else {
299: if (is->isperm) {
300: PetscViewerASCIISynchronizedPrintf(viewer,"Index set is permutation\n");
301: }
302: PetscViewerASCIISynchronizedPrintf(viewer,"Number of indices in set %D\n",n);
303: for (i=0; i<n; i++) {
304: PetscViewerASCIISynchronizedPrintf(viewer,"%D %D\n",i,idx[i]);
305: }
306: }
307: PetscViewerFlush(viewer);
308: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
309: } else if (isbinary) {
310: ISView_General_Binary(is,viewer);
311: }
312: return(0);
313: }
317: PetscErrorCode ISSort_General(IS is)
318: {
319: IS_General *sub = (IS_General*)is->data;
323: if (sub->sorted) return(0);
324: PetscSortInt(sub->n,sub->idx);
325: sub->sorted = PETSC_TRUE;
326: return(0);
327: }
331: PetscErrorCode ISSorted_General(IS is,PetscBool *flg)
332: {
333: IS_General *sub = (IS_General*)is->data;
336: *flg = sub->sorted;
337: return(0);
338: }
342: PetscErrorCode ISToGeneral_General(IS is)
343: {
345: return(0);
346: }
348: static struct _ISOps myops = { ISGetSize_General,
349: ISGetLocalSize_General,
350: ISGetIndices_General,
351: ISRestoreIndices_General,
352: ISInvertPermutation_General,
353: ISSort_General,
354: ISSorted_General,
355: ISDuplicate_General,
356: ISDestroy_General,
357: ISView_General,
358: ISIdentity_General,
359: ISCopy_General,
360: ISToGeneral_General,
361: ISOnComm_General,
362: ISSetBlockSize_General,
363: ISContiguousLocal_General};
367: PetscErrorCode ISCreateGeneral_Private(IS is)
368: {
370: IS_General *sub = (IS_General*)is->data;
371: PetscInt n = sub->n,i,min,max;
372: const PetscInt *idx = sub->idx;
373: PetscBool sorted = PETSC_TRUE;
374: PetscBool flg = PETSC_FALSE;
377: MPI_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)is));
378: for (i=1; i<n; i++) {
379: if (idx[i] < idx[i-1]) {sorted = PETSC_FALSE; break;}
380: }
381: if (n) min = max = idx[0];
382: else min = max = 0;
383: for (i=1; i<n; i++) {
384: if (idx[i] < min) min = idx[i];
385: if (idx[i] > max) max = idx[i];
386: }
387: sub->sorted = sorted;
388: is->min = min;
389: is->max = max;
390: is->isperm = PETSC_FALSE;
391: is->isidentity = PETSC_FALSE;
392: PetscOptionsGetBool(NULL,"-is_view",&flg,NULL);
393: if (flg) {
394: PetscViewer viewer;
395: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)is),&viewer);
396: ISView(is,viewer);
397: }
398: return(0);
399: }
403: /*@
404: ISCreateGeneral - Creates a data structure for an index set
405: containing a list of integers.
407: Collective on MPI_Comm
409: Input Parameters:
410: + comm - the MPI communicator
411: . n - the length of the index set
412: . idx - the list of integers
413: - mode - see PetscCopyMode for meaning of this flag.
415: Output Parameter:
416: . is - the new index set
418: Notes:
419: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
420: conceptually the same as MPI_Group operations. The IS are then
421: distributed sets of indices and thus certain operations on them are
422: collective.
425: Level: beginner
427: Concepts: index sets^creating
428: Concepts: IS^creating
430: .seealso: ISCreateStride(), ISCreateBlock(), ISAllGather()
431: @*/
432: PetscErrorCode ISCreateGeneral(MPI_Comm comm,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
433: {
437: ISCreate(comm,is);
438: ISSetType(*is,ISGENERAL);
439: ISGeneralSetIndices(*is,n,idx,mode);
440: return(0);
441: }
445: /*@
446: ISGeneralSetIndices - Sets the indices for an ISGENERAL index set
448: Collective on IS
450: Input Parameters:
451: + is - the index set
452: . n - the length of the index set
453: . idx - the list of integers
454: - mode - see PetscCopyMode for meaning of this flag.
456: Level: beginner
458: Concepts: index sets^creating
459: Concepts: IS^creating
461: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
462: @*/
463: PetscErrorCode ISGeneralSetIndices(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
464: {
468: PetscUseMethod(is,"ISGeneralSetIndices_C",(IS,PetscInt,const PetscInt[],PetscCopyMode),(is,n,idx,mode));
469: return(0);
470: }
474: PetscErrorCode ISGeneralSetIndices_General(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
475: {
477: IS_General *sub = (IS_General*)is->data;
480: if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
483: if (sub->allocated) {PetscFree(sub->idx);}
484: sub->n = n;
485: if (mode == PETSC_COPY_VALUES) {
486: PetscMalloc(n*sizeof(PetscInt),&sub->idx);
487: PetscLogObjectMemory(is,n*sizeof(PetscInt));
488: PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));
489: sub->allocated = PETSC_TRUE;
490: } else if (mode == PETSC_OWN_POINTER) {
491: sub->idx = (PetscInt*)idx;
492: sub->allocated = PETSC_TRUE;
493: } else {
494: sub->idx = (PetscInt*)idx;
495: sub->allocated = PETSC_FALSE;
496: }
497: ISCreateGeneral_Private(is);
498: return(0);
499: }
503: PETSC_EXTERN PetscErrorCode ISCreate_General(IS is)
504: {
506: IS_General *sub;
509: PetscMemcpy(is->ops,&myops,sizeof(myops));
510: PetscNewLog(is,IS_General,&sub);
511: is->data = (void*)sub;
512: is->bs = 1;
513: PetscObjectComposeFunction((PetscObject)is,"ISGeneralSetIndices_C",ISGeneralSetIndices_General);
514: return(0);
515: }