Actual source code: general.c

petsc-3.12.2 2019-11-22
Report Typos and Errors

  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>
  6:  #include <petsc/private/viewerimpl.h>
  7: #include <petsc/private/viewerhdf5impl.h>

  9: static PetscErrorCode ISDuplicate_General(IS is,IS *newIS)
 10: {
 12:   IS_General     *sub = (IS_General*)is->data;
 13:   PetscInt       n;

 16:   PetscLayoutGetLocalSize(is->map, &n);
 17:   ISCreateGeneral(PetscObjectComm((PetscObject) is), n, sub->idx, PETSC_COPY_VALUES, newIS);
 18:   return(0);
 19: }

 21: static PetscErrorCode ISDestroy_General(IS is)
 22: {
 23:   IS_General     *is_general = (IS_General*)is->data;

 27:   if (is_general->allocated) {PetscFree(is_general->idx);}
 28:   PetscObjectComposeFunction((PetscObject)is,"ISGeneralSetIndices_C",NULL);
 29:   PetscObjectComposeFunction((PetscObject)is,"ISGeneralFilter_C",NULL);
 30:   PetscFree(is->data);
 31:   return(0);
 32: }

 34: static PetscErrorCode ISIdentity_General(IS is, PetscBool *ident)
 35: {
 36:   IS_General *is_general = (IS_General*)is->data;
 37:   PetscInt   i,n,*idx = is_general->idx;

 41:   PetscLayoutGetLocalSize(is->map, &n);
 42:   is->isidentity = PETSC_TRUE;
 43:   *ident         = PETSC_TRUE;
 44:   for (i=0; i<n; i++) {
 45:     if (idx[i] != i) {
 46:       is->isidentity = PETSC_FALSE;
 47:       *ident         = PETSC_FALSE;
 48:       break;
 49:     }
 50:   }
 51:   return(0);
 52: }

 54: static PetscErrorCode ISCopy_General(IS is,IS isy)
 55: {
 56:   IS_General     *is_general = (IS_General*)is->data,*isy_general = (IS_General*)isy->data;
 57:   PetscInt       n, N, ny, Ny;

 61:   PetscLayoutGetLocalSize(is->map, &n);
 62:   PetscLayoutGetSize(is->map, &N);
 63:   PetscLayoutGetLocalSize(isy->map, &ny);
 64:   PetscLayoutGetSize(isy->map, &Ny);
 65:   if (n != ny || N != Ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Index sets incompatible");
 66:   isy_general->sorted = is_general->sorted;
 67:   PetscArraycpy(isy_general->idx,is_general->idx,n);
 68:   return(0);
 69: }

 71: static PetscErrorCode ISOnComm_General(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
 72: {
 74:   IS_General     *sub = (IS_General*)is->data;
 75:   PetscInt       n;

 78:   if (mode == PETSC_OWN_POINTER) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Cannot use PETSC_OWN_POINTER");
 79:   PetscLayoutGetLocalSize(is->map, &n);
 80:   ISCreateGeneral(comm,n,sub->idx,mode,newis);
 81:   return(0);
 82: }

 84: static PetscErrorCode ISSetBlockSize_General(IS is,PetscInt bs)
 85: {

 89:   PetscLayoutSetBlockSize(is->map, bs);
 90:   return(0);
 91: }

 93: static PetscErrorCode ISContiguousLocal_General(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
 94: {
 95:   IS_General *sub = (IS_General*)is->data;
 96:   PetscInt   n,i,p;

100:   *start  = 0;
101:   *contig = PETSC_TRUE;
102:   PetscLayoutGetLocalSize(is->map, &n);
103:   if (!n) return(0);
104:   p = sub->idx[0];
105:   if (p < gstart) goto nomatch;
106:   *start = p - gstart;
107:   if (n > gend-p) goto nomatch;
108:   for (i=1; i<n; i++,p++) {
109:     if (sub->idx[i] != p+1) goto nomatch;
110:   }
111:   return(0);
112: nomatch:
113:   *start  = -1;
114:   *contig = PETSC_FALSE;
115:   return(0);
116: }

118: static PetscErrorCode ISLocate_General(IS is,PetscInt key,PetscInt *location)
119: {
120:   IS_General     *sub = (IS_General*)is->data;
121:   PetscInt       numIdx, i;

125:   PetscLayoutGetLocalSize(is->map,&numIdx);
126:   if (sub->sorted) { PetscFindInt(key,numIdx,sub->idx,location);}
127:   else {
128:     const PetscInt *idx = sub->idx;

130:     *location = -1;
131:     for (i = 0; i < numIdx; i++) {
132:       if (idx[i] == key) {
133:         *location = i;
134:         return(0);
135:       }
136:     }
137:   }
138:   return(0);
139: }

141: static PetscErrorCode ISGetIndices_General(IS in,const PetscInt *idx[])
142: {
143:   IS_General *sub = (IS_General*)in->data;

146:   *idx = sub->idx;
147:   return(0);
148: }

150: static PetscErrorCode ISRestoreIndices_General(IS in,const PetscInt *idx[])
151: {
152:   IS_General *sub = (IS_General*)in->data;

155:   if (*idx != sub->idx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
156:   return(0);
157: }

159: static PetscErrorCode ISInvertPermutation_General(IS is,PetscInt nlocal,IS *isout)
160: {
161:   IS_General     *sub = (IS_General*)is->data;
162:   PetscInt       i,*ii,n,nstart;
163:   const PetscInt *idx = sub->idx;
164:   PetscMPIInt    size;
165:   IS             istmp,nistmp;

169:   PetscLayoutGetLocalSize(is->map, &n);
170:   MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
171:   if (size == 1) {
172:     PetscMalloc1(n,&ii);
173:     for (i=0; i<n; i++) ii[idx[i]] = i;
174:     ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_OWN_POINTER,isout);
175:     ISSetPermutation(*isout);
176:   } else {
177:     /* crude, nonscalable get entire IS on each processor */
178:     ISAllGather(is,&istmp);
179:     ISSetPermutation(istmp);
180:     ISInvertPermutation(istmp,PETSC_DECIDE,&nistmp);
181:     ISDestroy(&istmp);
182:     /* get the part we need */
183:     if (nlocal == PETSC_DECIDE) nlocal = n;
184:     MPI_Scan(&nlocal,&nstart,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)is));
185: #if defined(PETSC_USE_DEBUG)
186:     {
187:       PetscInt    N;
188:       PetscMPIInt rank;
189:       MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
190:       PetscLayoutGetSize(is->map, &N);
191:       if (rank == size-1) {
192:         if (nstart != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Sum of nlocal lengths %d != total IS length %d",nstart,N);
193:       }
194:     }
195: #endif
196:     nstart -= nlocal;
197:     ISGetIndices(nistmp,&idx);
198:     ISCreateGeneral(PetscObjectComm((PetscObject)is),nlocal,idx+nstart,PETSC_COPY_VALUES,isout);
199:     ISRestoreIndices(nistmp,&idx);
200:     ISDestroy(&nistmp);
201:   }
202:   return(0);
203: }

205: #if defined(PETSC_HAVE_HDF5)
206: static PetscErrorCode ISView_General_HDF5(IS is, PetscViewer viewer)
207: {
208:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
209:   hid_t           filespace;  /* file dataspace identifier */
210:   hid_t           chunkspace; /* chunk dataset property identifier */
211:   hid_t           dset_id;    /* dataset identifier */
212:   hid_t           memspace;   /* memory dataspace identifier */
213:   hid_t           inttype;    /* int type (H5T_NATIVE_INT or H5T_NATIVE_LLONG) */
214:   hid_t           file_id, group;
215:   hsize_t         dim, maxDims[3], dims[3], chunkDims[3], count[3],offset[3];
216:   PetscInt        bs, N, n, timestep, low;
217:   const PetscInt *ind;
218:   const char     *isname;
219:   PetscErrorCode  ierr;

222:   ISGetBlockSize(is,&bs);
223:   bs   = PetscMax(bs, 1); /* If N = 0, bs  = 0 as well */
224:   PetscViewerHDF5OpenGroup(viewer, &file_id, &group);
225:   PetscViewerHDF5GetTimestep(viewer, &timestep);

227:   /* Create the dataspace for the dataset.
228:    *
229:    * dims - holds the current dimensions of the dataset
230:    *
231:    * maxDims - holds the maximum dimensions of the dataset (unlimited
232:    * for the number of time steps with the current dimensions for the
233:    * other dimensions; so only additional time steps can be added).
234:    *
235:    * chunkDims - holds the size of a single time step (required to
236:    * permit extending dataset).
237:    */
238:   dim = 0;
239:   if (timestep >= 0) {
240:     dims[dim]      = timestep+1;
241:     maxDims[dim]   = H5S_UNLIMITED;
242:     chunkDims[dim] = 1;
243:     ++dim;
244:   }
245:   ISGetSize(is, &N);
246:   ISGetLocalSize(is, &n);
247:   PetscHDF5IntCast(N/bs,dims + dim);

249:   maxDims[dim]   = dims[dim];
250:   chunkDims[dim] = PetscMax(1,dims[dim]);
251:   ++dim;
252:   if (bs >= 1) {
253:     dims[dim]      = bs;
254:     maxDims[dim]   = dims[dim];
255:     chunkDims[dim] = dims[dim];
256:     ++dim;
257:   }
258:   PetscStackCallHDF5Return(filespace,H5Screate_simple,(dim, dims, maxDims));

260: #if defined(PETSC_USE_64BIT_INDICES)
261:   inttype = H5T_NATIVE_LLONG;
262: #else
263:   inttype = H5T_NATIVE_INT;
264: #endif

266:   /* Create the dataset with default properties and close filespace */
267:   PetscObjectGetName((PetscObject) is, &isname);
268:   if (!H5Lexists(group, isname, H5P_DEFAULT)) {
269:     /* Create chunk */
270:     PetscStackCallHDF5Return(chunkspace,H5Pcreate,(H5P_DATASET_CREATE));
271:     PetscStackCallHDF5(H5Pset_chunk,(chunkspace, dim, chunkDims));

273:     PetscStackCallHDF5Return(dset_id,H5Dcreate2,(group, isname, inttype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
274:     PetscStackCallHDF5(H5Pclose,(chunkspace));
275:   } else {
276:     PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, isname, H5P_DEFAULT));
277:     PetscStackCallHDF5(H5Dset_extent,(dset_id, dims));
278:   }
279:   PetscStackCallHDF5(H5Sclose,(filespace));

281:   /* Each process defines a dataset and writes it to the hyperslab in the file */
282:   dim = 0;
283:   if (timestep >= 0) {
284:     count[dim] = 1;
285:     ++dim;
286:   }
287:   PetscHDF5IntCast(n/bs,count + dim);
288:   ++dim;
289:   if (bs >= 1) {
290:     count[dim] = bs;
291:     ++dim;
292:   }
293:   if (n > 0) {
294:     PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL));
295:   } else {
296:     /* Can't create dataspace with zero for any dimension, so create null dataspace. */
297:     PetscStackCallHDF5Return(memspace,H5Screate,(H5S_NULL));
298:   }

300:   /* Select hyperslab in the file */
301:   PetscLayoutGetRange(is->map, &low, NULL);
302:   dim  = 0;
303:   if (timestep >= 0) {
304:     offset[dim] = timestep;
305:     ++dim;
306:   }
307:   PetscHDF5IntCast(low/bs,offset + dim);
308:   ++dim;
309:   if (bs >= 1) {
310:     offset[dim] = 0;
311:     ++dim;
312:   }
313:   if (n > 0) {
314:     PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id));
315:     PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
316:   } else {
317:     /* Create null filespace to match null memspace. */
318:     PetscStackCallHDF5Return(filespace,H5Screate,(H5S_NULL));
319:   }

321:   ISGetIndices(is, &ind);
322:   PetscStackCallHDF5(H5Dwrite,(dset_id, inttype, memspace, filespace, hdf5->dxpl_id, ind));
323:   PetscStackCallHDF5(H5Fflush,(file_id, H5F_SCOPE_GLOBAL));
324:   ISRestoreIndices(is, &ind);

326:   /* Close/release resources */
327:   PetscStackCallHDF5(H5Gclose,(group));
328:   PetscStackCallHDF5(H5Sclose,(filespace));
329:   PetscStackCallHDF5(H5Sclose,(memspace));
330:   PetscStackCallHDF5(H5Dclose,(dset_id));
331:   PetscInfo1(is, "Wrote IS object with name %s\n", isname);
332:   return(0);
333: }
334: #endif

336: static PetscErrorCode ISView_General_Binary(IS is,PetscViewer viewer)
337: {
339:   PetscBool      skipHeader,useMPIIO;
340:   IS_General     *isa = (IS_General*) is->data;
341:   PetscMPIInt    rank,size,mesgsize,tag = ((PetscObject)viewer)->tag, mesglen;
342:   PetscInt       n,N,len,j,tr[2];
343:   int            fdes;
344:   MPI_Status     status;
345:   PetscInt       message_count,flowcontrolcount,*values;

348:   /* ISGetLayout(is,&map); */
349:   PetscLayoutGetLocalSize(is->map, &n);
350:   PetscLayoutGetSize(is->map, &N);

352:   tr[0] = IS_FILE_CLASSID;
353:   tr[1] = N;

355:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
356:   if (!skipHeader) {
357:     PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_FALSE);
358:   }

360:   PetscViewerBinaryGetUseMPIIO(viewer,&useMPIIO);
361: #if defined(PETSC_HAVE_MPIIO)
362:   if (useMPIIO) {
363:     MPI_File       mfdes;
364:     MPI_Offset     off;
365:     PetscMPIInt    lsize;
366:     PetscInt       rstart;
367:     const PetscInt *iarray;

369:     PetscMPIIntCast(n,&lsize);
370:     PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);
371:     PetscViewerBinaryGetMPIIOOffset(viewer,&off);
372:     PetscLayoutGetRange(is->map,&rstart,NULL);
373:     off += rstart*(MPI_Offset)sizeof(PetscInt); /* off is MPI_Offset, not PetscMPIInt */
374:     MPI_File_set_view(mfdes,off,MPIU_INT,MPIU_INT,(char*)"native",MPI_INFO_NULL);
375:     ISGetIndices(is,&iarray);
376:     MPIU_File_write_all(mfdes,(void*)iarray,lsize,MPIU_INT,MPI_STATUS_IGNORE);
377:     ISRestoreIndices(is,&iarray);
378:     PetscViewerBinaryAddMPIIOOffset(viewer,N*(MPI_Offset)sizeof(PetscInt));
379:     return(0);
380:   }
381: #endif

383:   PetscViewerBinaryGetDescriptor(viewer,&fdes);
384:   MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
385:   MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);

387:   /* determine maximum message to arrive */
388:   MPI_Reduce(&n,&len,1,MPIU_INT,MPI_MAX,0,PetscObjectComm((PetscObject)is));

390:   PetscViewerFlowControlStart(viewer,&message_count,&flowcontrolcount);
391:   if (!rank) {
392:     PetscBinaryWrite(fdes,isa->idx,n,PETSC_INT,PETSC_FALSE);

394:     PetscMalloc1(len,&values);
395:     PetscMPIIntCast(len,&mesgsize);
396:     /* receive and save messages */
397:     for (j=1; j<size; j++) {
398:       PetscViewerFlowControlStepMaster(viewer,j,&message_count,flowcontrolcount);
399:       MPI_Recv(values,mesgsize,MPIU_INT,j,tag,PetscObjectComm((PetscObject)is),&status);
400:       MPI_Get_count(&status,MPIU_INT,&mesglen);
401:       PetscBinaryWrite(fdes,values,(PetscInt)mesglen,PETSC_INT,PETSC_TRUE);
402:     }
403:     PetscViewerFlowControlEndMaster(viewer,&message_count);
404:     PetscFree(values);
405:   } else {
406:     PetscViewerFlowControlStepWorker(viewer,rank,&message_count);
407:     PetscMPIIntCast(n,&mesgsize);
408:     MPI_Send(isa->idx,mesgsize,MPIU_INT,0,tag,PetscObjectComm((PetscObject)is));
409:     PetscViewerFlowControlEndWorker(viewer,&message_count);
410:   }
411:   return(0);
412: }

414: static PetscErrorCode ISView_General(IS is,PetscViewer viewer)
415: {
416:   IS_General     *sub = (IS_General*)is->data;
418:   PetscInt       i,n,*idx = sub->idx;
419:   PetscBool      iascii,isbinary,ishdf5;

422:   PetscLayoutGetLocalSize(is->map, &n);
423:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
424:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
425:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
426:   if (iascii) {
427:     MPI_Comm          comm;
428:     PetscMPIInt       rank,size;
429:     PetscViewerFormat fmt;

431:     PetscObjectGetComm((PetscObject)viewer,&comm);
432:     MPI_Comm_rank(comm,&rank);
433:     MPI_Comm_size(comm,&size);

435:     PetscViewerGetFormat(viewer,&fmt);
436:     PetscViewerASCIIPushSynchronized(viewer);
437:     if (size > 1) {
438:       if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
439:         const char* name;

441:         PetscObjectGetName((PetscObject)is,&name);
442:         PetscViewerASCIISynchronizedPrintf(viewer,"%s_%d = [...\n",name,rank);
443:         for (i=0; i<n; i++) {
444:           PetscViewerASCIISynchronizedPrintf(viewer,"%D\n",idx[i]+1);
445:         }
446:         PetscViewerASCIISynchronizedPrintf(viewer,"];\n");
447:       } else {
448:         PetscInt st = 0;

450:         if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
451:         if (is->isperm) {
452:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutation\n",rank);
453:         }
454:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in set %D\n",rank,n);
455:         for (i=0; i<n; i++) {
456:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i + st,idx[i]);
457:         }
458:       }
459:     } else {
460:       if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
461:         const char* name;

463:         PetscObjectGetName((PetscObject)is,&name);
464:         PetscViewerASCIISynchronizedPrintf(viewer,"%s = [...\n",name);
465:         for (i=0; i<n; i++) {
466:           PetscViewerASCIISynchronizedPrintf(viewer,"%D\n",idx[i]+1);
467:         }
468:         PetscViewerASCIISynchronizedPrintf(viewer,"];\n");
469:       } else {
470:         PetscInt st = 0;

472:         if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
473:         if (is->isperm) {
474:           PetscViewerASCIISynchronizedPrintf(viewer,"Index set is permutation\n");
475:         }
476:         PetscViewerASCIISynchronizedPrintf(viewer,"Number of indices in set %D\n",n);
477:         for (i=0; i<n; i++) {
478:           PetscViewerASCIISynchronizedPrintf(viewer,"%D %D\n",i + st,idx[i]);
479:         }
480:       }
481:     }
482:     PetscViewerFlush(viewer);
483:     PetscViewerASCIIPopSynchronized(viewer);
484:   } else if (isbinary) {
485:     ISView_General_Binary(is,viewer);
486:   } else if (ishdf5) {
487: #if defined(PETSC_HAVE_HDF5)
488:     ISView_General_HDF5(is,viewer);
489: #endif
490:   }
491:   return(0);
492: }

494: static PetscErrorCode ISSort_General(IS is)
495: {
496:   IS_General     *sub = (IS_General*)is->data;
497:   PetscInt       n;

501:   if (sub->sorted) return(0);
502:   PetscLayoutGetLocalSize(is->map, &n);
503:   PetscSortInt(n,sub->idx);
504:   sub->sorted = PETSC_TRUE;
505:   return(0);
506: }

508: static PetscErrorCode ISSortRemoveDups_General(IS is)
509: {
510:   IS_General     *sub = (IS_General*)is->data;
511:   PetscLayout    map;
512:   PetscInt       n;

516:   PetscLayoutGetLocalSize(is->map, &n);
517:   if (sub->sorted) {
518:     PetscSortedRemoveDupsInt(&n,sub->idx);
519:   } else {
520:     PetscSortRemoveDupsInt(&n,sub->idx);
521:   }
522:   PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map);
523:   PetscLayoutDestroy(&is->map);
524:   is->map = map;
525:   sub->sorted = PETSC_TRUE;
526:   return(0);
527: }

529: static PetscErrorCode ISSorted_General(IS is,PetscBool  *flg)
530: {
531:   IS_General *sub = (IS_General*)is->data;

534:   *flg = sub->sorted;
535:   return(0);
536: }

538: PetscErrorCode  ISToGeneral_General(IS is)
539: {
541:   return(0);
542: }

544: static struct _ISOps myops = { ISGetIndices_General,
545:                                ISRestoreIndices_General,
546:                                ISInvertPermutation_General,
547:                                ISSort_General,
548:                                ISSortRemoveDups_General,
549:                                ISSorted_General,
550:                                ISDuplicate_General,
551:                                ISDestroy_General,
552:                                ISView_General,
553:                                ISLoad_Default,
554:                                ISIdentity_General,
555:                                ISCopy_General,
556:                                ISToGeneral_General,
557:                                ISOnComm_General,
558:                                ISSetBlockSize_General,
559:                                ISContiguousLocal_General,
560:                                ISLocate_General};

562: PETSC_INTERN PetscErrorCode ISSetUp_General(IS);

564: PetscErrorCode ISSetUp_General(IS is)
565: {
567:   IS_General     *sub = (IS_General*)is->data;
568:   const PetscInt *idx = sub->idx;
569:   PetscInt       n,i,min,max;

572:   PetscLayoutGetLocalSize(is->map, &n);

574:   sub->sorted = PETSC_TRUE;
575:   for (i=1; i<n; i++) {
576:     if (idx[i] < idx[i-1]) {sub->sorted = PETSC_FALSE; break;}
577:   }
578:   if (n) {
579:     min = max = idx[0];
580:     for (i=1; i<n; i++) {
581:       if (idx[i] < min) min = idx[i];
582:       if (idx[i] > max) max = idx[i];
583:     }
584:     is->min = min;
585:     is->max = max;
586:   } else {
587:     is->min = PETSC_MAX_INT;
588:     is->max = PETSC_MIN_INT;
589:   }
590:   is->isperm     = PETSC_FALSE;
591:   is->isidentity = PETSC_FALSE;
592:   return(0);
593: }

595: /*@
596:    ISCreateGeneral - Creates a data structure for an index set
597:    containing a list of integers.

599:    Collective

601:    Input Parameters:
602: +  comm - the MPI communicator
603: .  n - the length of the index set
604: .  idx - the list of integers
605: -  mode - PETSC_COPY_VALUES, PETSC_OWN_POINTER, or PETSC_USE_POINTER; see PetscCopyMode for meaning of this flag.

607:    Output Parameter:
608: .  is - the new index set

610:    Notes:
611:    When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
612:    conceptually the same as MPI_Group operations. The IS are then
613:    distributed sets of indices and thus certain operations on them are
614:    collective.


617:    Level: beginner


620: .seealso: ISCreateStride(), ISCreateBlock(), ISAllGather(), PETSC_COPY_VALUES, PETSC_OWN_POINTER, PETSC_USE_POINTER, PetscCopyMode
621: @*/
622: PetscErrorCode  ISCreateGeneral(MPI_Comm comm,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
623: {

627:   ISCreate(comm,is);
628:   ISSetType(*is,ISGENERAL);
629:   ISGeneralSetIndices(*is,n,idx,mode);
630:   return(0);
631: }

633: /*@
634:    ISGeneralSetIndices - Sets the indices for an ISGENERAL index set

636:    Collective on IS

638:    Input Parameters:
639: +  is - the index set
640: .  n - the length of the index set
641: .  idx - the list of integers
642: -  mode - see PetscCopyMode for meaning of this flag.

644:    Level: beginner


647: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
648: @*/
649: PetscErrorCode  ISGeneralSetIndices(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
650: {

654:   PetscUseMethod(is,"ISGeneralSetIndices_C",(IS,PetscInt,const PetscInt[],PetscCopyMode),(is,n,idx,mode));
655:   return(0);
656: }

658: PetscErrorCode  ISGeneralSetIndices_General(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
659: {
660:   PetscLayout    map;
662:   IS_General     *sub = (IS_General*)is->data;

665:   if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");

668:   PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is),n,PETSC_DECIDE,is->map->bs,&map);
669:   PetscLayoutDestroy(&is->map);
670:   is->map = map;

672:   if (sub->allocated) {PetscFree(sub->idx);}
673:   if (mode == PETSC_COPY_VALUES) {
674:     PetscMalloc1(n,&sub->idx);
675:     PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
676:     PetscArraycpy(sub->idx,idx,n);
677:     sub->allocated = PETSC_TRUE;
678:   } else if (mode == PETSC_OWN_POINTER) {
679:     sub->idx = (PetscInt*)idx;
680:     PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
681:     sub->allocated = PETSC_TRUE;
682:   } else {
683:     sub->idx = (PetscInt*)idx;
684:     sub->allocated = PETSC_FALSE;
685:   }

687:   ISSetUp_General(is);
688:   ISViewFromOptions(is,NULL,"-is_view");
689:   return(0);
690: }

692: static PetscErrorCode ISGeneralFilter_General(IS is, PetscInt start, PetscInt end)
693: {
694:   IS_General     *sub = (IS_General*)is->data;
695:   PetscInt       *idx = sub->idx,*idxnew;
696:   PetscInt       i,n = is->map->n,nnew = 0,o;

700:   for (i=0; i<n; ++i)
701:     if (idx[i] >= start && idx[i] < end)
702:       nnew++;
703:   PetscMalloc1(nnew, &idxnew);
704:   for (o=0, i=0; i<n; i++) {
705:     if (idx[i] >= start && idx[i] < end)
706:       idxnew[o++] = idx[i];
707:   }
708:   ISGeneralSetIndices_General(is,nnew,idxnew,PETSC_OWN_POINTER);
709:   return(0);
710: }

712: /*@
713:    ISGeneralFilter - Remove all points outside of [start, end)

715:    Collective on IS

717:    Input Parameters:
718: +  is - the index set
719: .  start - the lowest index kept
720: -  end - one more than the highest index kept

722:    Level: beginner

724: .seealso: ISCreateGeneral(), ISGeneralSetIndices()
725: @*/
726: PetscErrorCode ISGeneralFilter(IS is, PetscInt start, PetscInt end)
727: {

732:   PetscUseMethod(is,"ISGeneralFilter_C",(IS,PetscInt,PetscInt),(is,start,end));
733:   return(0);
734: }

736: PETSC_EXTERN PetscErrorCode ISCreate_General(IS is)
737: {
739:   IS_General     *sub;

742:   PetscNewLog(is,&sub);
743:   is->data = (void *) sub;
744:   PetscMemcpy(is->ops,&myops,sizeof(myops));
745:   PetscObjectComposeFunction((PetscObject)is,"ISGeneralSetIndices_C",ISGeneralSetIndices_General);
746:   PetscObjectComposeFunction((PetscObject)is,"ISGeneralFilter_C",ISGeneralFilter_General);
747:   return(0);
748: }