Actual source code: general.c

petsc-3.4.4 2014-03-13
  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: }