Actual source code: general.c
petsc-3.3-p7 2013-05-11
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> /*I "petscis.h" I*/
6: #include <petscvec.h>
10: PetscErrorCode ISDuplicate_General(IS is,IS *newIS)
11: {
13: IS_General *sub = (IS_General *)is->data;
16: ISCreateGeneral(((PetscObject)is)->comm,sub->n,sub->idx,PETSC_COPY_VALUES,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: PetscObjectComposeFunctionDynamic((PetscObject)is,"ISGeneralSetIndices_C","",0);
32: PetscFree(is->data);
33: return(0);
34: }
38: PetscErrorCode ISIdentity_General(IS is,PetscBool *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: static PetscErrorCode ISCopy_General(IS is,IS isy)
59: {
60: IS_General *is_general = (IS_General*)is->data,*isy_general = (IS_General*)isy->data;
64: if (is_general->n != isy_general->n || is_general->N != isy_general->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Index sets incompatible");
65: isy_general->sorted = is_general->sorted;
66: PetscMemcpy(isy_general->idx,is_general->idx,is_general->n*sizeof(PetscInt));
67: return(0);
68: }
72: PetscErrorCode ISOnComm_General(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
73: {
75: IS_General *sub = (IS_General*)is->data;
78: if (mode == PETSC_OWN_POINTER) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Cannot use PETSC_OWN_POINTER");
79: ISCreateGeneral(comm,sub->n,sub->idx,mode,newis);
80: return(0);
81: }
85: static PetscErrorCode ISSetBlockSize_General(IS is,PetscInt bs)
86: {
87: IS_General *sub = (IS_General*)is->data;
90: if (sub->N % bs) SETERRQ2(((PetscObject)is)->comm,PETSC_ERR_ARG_SIZ,"Block size %D does not divide global size %D",bs,sub->N);
91: if (sub->n % bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Block size %D does not divide local size %D",bs,sub->n);
92: #if defined(PETSC_USE_DEBUG)
93: {
94: PetscInt i,j;
95: for (i=0; i<sub->n; i+=bs) {
96: for (j=0; j<bs; j++) {
97: 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);
98: }
99: }
100: }
101: #endif
102: is->bs = bs;
103: return(0);
104: }
108: static PetscErrorCode ISContiguousLocal_General(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
109: {
110: IS_General *sub = (IS_General*)is->data;
111: PetscInt i,p;
114: *start = 0;
115: *contig = PETSC_TRUE;
116: if (!sub->n) return(0);
117: p = sub->idx[0];
118: if (p < gstart) goto nomatch;
119: *start = p - gstart;
120: if (sub->n > gend-p) goto nomatch;
121: for (i=1; i<sub->n; i++,p++) {
122: if (sub->idx[i] != p+1) goto nomatch;
123: }
124: return(0);
125: nomatch:
126: *start = -1;
127: *contig = PETSC_FALSE;
128: return(0);
129: }
133: PetscErrorCode ISGetIndices_General(IS in,const PetscInt *idx[])
134: {
135: IS_General *sub = (IS_General*)in->data;
138: *idx = sub->idx;
139: return(0);
140: }
144: PetscErrorCode ISRestoreIndices_General(IS in,const PetscInt *idx[])
145: {
146: IS_General *sub = (IS_General*)in->data;
149: if (*idx != sub->idx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
150: return(0);
151: }
155: PetscErrorCode ISGetSize_General(IS is,PetscInt *size)
156: {
157: IS_General *sub = (IS_General *)is->data;
160: *size = sub->N;
161: return(0);
162: }
166: PetscErrorCode ISGetLocalSize_General(IS is,PetscInt *size)
167: {
168: IS_General *sub = (IS_General *)is->data;
171: *size = sub->n;
172: return(0);
173: }
177: PetscErrorCode ISInvertPermutation_General(IS is,PetscInt nlocal,IS *isout)
178: {
179: IS_General *sub = (IS_General *)is->data;
180: PetscInt i,*ii,n = sub->n,nstart;
181: const PetscInt *idx = sub->idx;
182: PetscMPIInt size;
183: IS istmp,nistmp;
187: MPI_Comm_size(((PetscObject)is)->comm,&size);
188: if (size == 1) {
189: PetscMalloc(n*sizeof(PetscInt),&ii);
190: for (i=0; i<n; i++) {
191: ii[idx[i]] = i;
192: }
193: ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_OWN_POINTER,isout);
194: ISSetPermutation(*isout);
195: } else {
196: /* crude, nonscalable get entire IS on each processor */
197: if (nlocal == PETSC_DECIDE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Do not yet support nlocal of PETSC_DECIDE");
198: ISAllGather(is,&istmp);
199: ISSetPermutation(istmp);
200: ISInvertPermutation(istmp,PETSC_DECIDE,&nistmp);
201: ISDestroy(&istmp);
202: /* get the part we need */
203: MPI_Scan(&nlocal,&nstart,1,MPIU_INT,MPI_SUM,((PetscObject)is)->comm);
204: #if defined(PETSC_USE_DEBUG)
205: {
206: PetscMPIInt rank;
207: MPI_Comm_rank(((PetscObject)is)->comm,&rank);
208: if (rank == size-1) {
209: if (nstart != sub->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Sum of nlocal lengths %d != total IS length %d",nstart,sub->N);
210: }
211: }
212: #endif
213: nstart -= nlocal;
214: ISGetIndices(nistmp,&idx);
215: ISCreateGeneral(((PetscObject)is)->comm,nlocal,idx+nstart,PETSC_COPY_VALUES,isout);
216: ISRestoreIndices(nistmp,&idx);
217: ISDestroy(&nistmp);
218: }
219: return(0);
220: }
224: PetscErrorCode ISView_General_Binary(IS is,PetscViewer viewer)
225: {
227: IS_General *isa = (IS_General*) is->data;
228: PetscMPIInt rank,size,mesgsize,tag = ((PetscObject)viewer)->tag, mesglen;
229: PetscInt len,j,tr[2];
230: int fdes;
231: MPI_Status status;
232: PetscInt message_count,flowcontrolcount,*values;
235: PetscViewerBinaryGetDescriptor(viewer,&fdes);
237: /* determine maximum message to arrive */
238: MPI_Comm_rank(((PetscObject)is)->comm,&rank);
239: MPI_Comm_size(((PetscObject)is)->comm,&size);
241: tr[0] = IS_FILE_CLASSID;
242: tr[1] = isa->N;
243: PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_FALSE);
244: MPI_Reduce(&isa->n,&len,1,MPIU_INT,MPI_SUM,0,((PetscObject)is)->comm);
246: PetscViewerFlowControlStart(viewer,&message_count,&flowcontrolcount);
247: if (!rank) {
248: PetscBinaryWrite(fdes,isa->idx,isa->n,PETSC_INT,PETSC_FALSE);
250: PetscMalloc(len*sizeof(PetscInt),&values);
251: mesgsize = PetscMPIIntCast(len);
252: /* receive and save messages */
253: for (j=1; j<size; j++) {
254: PetscViewerFlowControlStepMaster(viewer,j,message_count,flowcontrolcount);
255: MPI_Recv(values,mesgsize,MPIU_INT,j,tag,((PetscObject)is)->comm,&status);
256: MPI_Get_count(&status,MPIU_INT,&mesglen);
257: PetscBinaryWrite(fdes,values,(PetscInt)mesglen,PETSC_INT,PETSC_TRUE);
258: }
259: PetscViewerFlowControlEndMaster(viewer,message_count);
260: PetscFree(values);
261: } else {
262: PetscViewerFlowControlStepWorker(viewer,rank,message_count);
263: mesgsize = PetscMPIIntCast(isa->n);
264: MPI_Send(isa->idx,mesgsize,MPIU_INT,0,tag,((PetscObject)is)->comm);
265: PetscViewerFlowControlEndWorker(viewer,message_count);
266: }
267: return(0);
268: }
272: PetscErrorCode ISView_General(IS is,PetscViewer viewer)
273: {
274: IS_General *sub = (IS_General *)is->data;
276: PetscInt i,n = sub->n,*idx = sub->idx;
277: PetscBool iascii,isbinary;
280: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
281: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
282: if (iascii) {
283: MPI_Comm comm;
284: PetscMPIInt rank,size;
286: PetscObjectGetComm((PetscObject)viewer,&comm);
287: MPI_Comm_rank(comm,&rank);
288: MPI_Comm_size(comm,&size);
290: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
291: if (size > 1) {
292: if (is->isperm) {
293: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutation\n",rank);
294: }
295: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in set %D\n",rank,n);
296: for (i=0; i<n; i++) {
297: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,idx[i]);
298: }
299: } else {
300: if (is->isperm) {
301: PetscViewerASCIISynchronizedPrintf(viewer,"Index set is permutation\n");
302: }
303: PetscViewerASCIISynchronizedPrintf(viewer,"Number of indices in set %D\n",n);
304: for (i=0; i<n; i++) {
305: PetscViewerASCIISynchronizedPrintf(viewer,"%D %D\n",i,idx[i]);
306: }
307: }
308: PetscViewerFlush(viewer);
309: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
310: } else if (isbinary) {
311: ISView_General_Binary(is,viewer);
312: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for this object",((PetscObject)viewer)->type_name);
313: return(0);
314: }
318: PetscErrorCode ISSort_General(IS is)
319: {
320: IS_General *sub = (IS_General *)is->data;
324: if (sub->sorted) return(0);
325: PetscSortInt(sub->n,sub->idx);
326: sub->sorted = PETSC_TRUE;
327: return(0);
328: }
332: PetscErrorCode ISSorted_General(IS is,PetscBool *flg)
333: {
334: IS_General *sub = (IS_General *)is->data;
337: *flg = sub->sorted;
338: return(0);
339: }
343: PetscErrorCode ISToGeneral_General(IS is)
344: {
346: return(0);
347: }
349: static struct _ISOps myops = { ISGetSize_General,
350: ISGetLocalSize_General,
351: ISGetIndices_General,
352: ISRestoreIndices_General,
353: ISInvertPermutation_General,
354: ISSort_General,
355: ISSorted_General,
356: ISDuplicate_General,
357: ISDestroy_General,
358: ISView_General,
359: ISIdentity_General,
360: ISCopy_General,
361: ISToGeneral_General,
362: ISOnComm_General,
363: ISSetBlockSize_General,
364: ISContiguousLocal_General
365: };
369: PetscErrorCode ISCreateGeneral_Private(IS is)
370: {
372: IS_General *sub = (IS_General*)is->data;
373: PetscInt n = sub->n,i,min,max;
374: const PetscInt *idx = sub->idx;
375: PetscBool sorted = PETSC_TRUE;
376: PetscBool flg = PETSC_FALSE;
379: MPI_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,((PetscObject)is)->comm);
380: for (i=1; i<n; i++) {
381: if (idx[i] < idx[i-1]) {sorted = PETSC_FALSE; break;}
382: }
383: if (n) {min = max = idx[0];} else {min = max = 0;}
384: for (i=1; i<n; i++) {
385: if (idx[i] < min) min = idx[i];
386: if (idx[i] > max) max = idx[i];
387: }
388: sub->sorted = sorted;
389: is->min = min;
390: is->max = max;
391: is->isperm = PETSC_FALSE;
392: is->isidentity = PETSC_FALSE;
393: PetscOptionsGetBool(PETSC_NULL,"-is_view",&flg,PETSC_NULL);
394: if (flg) {
395: PetscViewer viewer;
396: PetscViewerASCIIGetStdout(((PetscObject)is)->comm,&viewer);
397: ISView(is,viewer);
398: }
399: return(0);
400: }
404: /*@
405: ISCreateGeneral - Creates a data structure for an index set
406: containing a list of integers.
408: Collective on MPI_Comm
410: Input Parameters:
411: + comm - the MPI communicator
412: . n - the length of the index set
413: . idx - the list of integers
414: - mode - see PetscCopyMode for meaning of this flag.
416: Output Parameter:
417: . is - the new index set
419: Notes:
420: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
421: conceptually the same as MPI_Group operations. The IS are then
422: distributed sets of indices and thus certain operations on them are
423: collective.
425:
426: Level: beginner
428: Concepts: index sets^creating
429: Concepts: IS^creating
431: .seealso: ISCreateStride(), ISCreateBlock(), ISAllGather()
432: @*/
433: PetscErrorCode ISCreateGeneral(MPI_Comm comm,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
434: {
438: ISCreate(comm,is);
439: ISSetType(*is,ISGENERAL);
440: ISGeneralSetIndices(*is,n,idx,mode);
441: return(0);
442: }
446: /*@
447: ISGeneralSetIndices - Sets the indices for an ISGENERAL index set
449: Collective on IS
451: Input Parameters:
452: + is - the index set
453: . n - the length of the index set
454: . idx - the list of integers
455: - mode - see PetscCopyMode for meaning of this flag.
457: Level: beginner
459: Concepts: index sets^creating
460: Concepts: IS^creating
462: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
463: @*/
464: PetscErrorCode ISGeneralSetIndices(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
465: {
468: PetscUseMethod(is,"ISGeneralSetIndices_C",(IS,PetscInt,const PetscInt[],PetscCopyMode),(is,n,idx,mode));
469: return(0);
470: }
472: EXTERN_C_BEGIN
475: PetscErrorCode ISGeneralSetIndices_General(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
476: {
478: IS_General *sub = (IS_General*)is->data;
481: if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
484: if (sub->allocated) {PetscFree(sub->idx);}
485: sub->n = n;
486: if (mode == PETSC_COPY_VALUES) {
487: PetscMalloc(n*sizeof(PetscInt),&sub->idx);
488: PetscLogObjectMemory(is,n*sizeof(PetscInt));
489: PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));
490: sub->allocated = PETSC_TRUE;
491: } else if (mode == PETSC_OWN_POINTER) {
492: sub->idx = (PetscInt*)idx;
493: sub->allocated = PETSC_TRUE;
494: } else {
495: sub->idx = (PetscInt*)idx;
496: sub->allocated = PETSC_FALSE;
497: }
498: ISCreateGeneral_Private(is);
499: return(0);
500: }
501: EXTERN_C_END
506: EXTERN_C_BEGIN
509: PetscErrorCode ISCreate_General(IS is)
510: {
512: IS_General *sub;
515: PetscMemcpy(is->ops,&myops,sizeof(myops));
516: PetscNewLog(is,IS_General,&sub);
517: is->data = (void*)sub;
518: is->bs = 1;
519: PetscObjectComposeFunctionDynamic((PetscObject)is,"ISGeneralSetIndices_C","ISGeneralSetIndices_General",ISGeneralSetIndices_General);
520: return(0);
521: }
522: EXTERN_C_END