Actual source code: general.c


  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/viewerhdf5impl.h>

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

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

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

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

 33: static PetscErrorCode ISCopy_General(IS is,IS isy)
 34: {
 35:   IS_General     *is_general = (IS_General*)is->data,*isy_general = (IS_General*)isy->data;
 36:   PetscInt       n, N, ny, Ny;

 40:   PetscLayoutGetLocalSize(is->map, &n);
 41:   PetscLayoutGetSize(is->map, &N);
 42:   PetscLayoutGetLocalSize(isy->map, &ny);
 43:   PetscLayoutGetSize(isy->map, &Ny);
 44:   if (n != ny || N != Ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Index sets incompatible");
 45:   PetscArraycpy(isy_general->idx,is_general->idx,n);
 46:   return(0);
 47: }

 49: static PetscErrorCode ISOnComm_General(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
 50: {
 52:   IS_General     *sub = (IS_General*)is->data;
 53:   PetscInt       n;

 56:   if (mode == PETSC_OWN_POINTER) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Cannot use PETSC_OWN_POINTER");
 57:   PetscLayoutGetLocalSize(is->map, &n);
 58:   ISCreateGeneral(comm,n,sub->idx,mode,newis);
 59:   return(0);
 60: }

 62: static PetscErrorCode ISSetBlockSize_General(IS is,PetscInt bs)
 63: {

 67:   PetscLayoutSetBlockSize(is->map, bs);
 68:   return(0);
 69: }

 71: static PetscErrorCode ISContiguousLocal_General(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
 72: {
 73:   IS_General *sub = (IS_General*)is->data;
 74:   PetscInt   n,i,p;

 78:   *start  = 0;
 79:   *contig = PETSC_TRUE;
 80:   PetscLayoutGetLocalSize(is->map, &n);
 81:   if (!n) return(0);
 82:   p = sub->idx[0];
 83:   if (p < gstart) goto nomatch;
 84:   *start = p - gstart;
 85:   if (n > gend-p) goto nomatch;
 86:   for (i=1; i<n; i++,p++) {
 87:     if (sub->idx[i] != p+1) goto nomatch;
 88:   }
 89:   return(0);
 90: nomatch:
 91:   *start  = -1;
 92:   *contig = PETSC_FALSE;
 93:   return(0);
 94: }

 96: static PetscErrorCode ISLocate_General(IS is,PetscInt key,PetscInt *location)
 97: {
 98:   IS_General     *sub = (IS_General*)is->data;
 99:   PetscInt       numIdx, i;
100:   PetscBool      sorted;

104:   PetscLayoutGetLocalSize(is->map,&numIdx);
105:   ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sorted);
106:   if (sorted) {PetscFindInt(key,numIdx,sub->idx,location);}
107:   else {
108:     const PetscInt *idx = sub->idx;

110:     *location = -1;
111:     for (i = 0; i < numIdx; i++) {
112:       if (idx[i] == key) {
113:         *location = i;
114:         return(0);
115:       }
116:     }
117:   }
118:   return(0);
119: }

121: static PetscErrorCode ISGetIndices_General(IS in,const PetscInt *idx[])
122: {
123:   IS_General *sub = (IS_General*)in->data;

126:   *idx = sub->idx;
127:   return(0);
128: }

130: static PetscErrorCode ISRestoreIndices_General(IS in,const PetscInt *idx[])
131: {
132:   IS_General *sub = (IS_General*)in->data;

135:    /* F90Array1dCreate() inside ISRestoreArrayF90() does not keep array when zero length array */
136:   if (in->map->n > 0  && *idx != sub->idx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
137:   return(0);
138: }

140: static PetscErrorCode ISInvertPermutation_General(IS is,PetscInt nlocal,IS *isout)
141: {
142:   IS_General     *sub = (IS_General*)is->data;
143:   PetscInt       i,*ii,n,nstart;
144:   const PetscInt *idx = sub->idx;
145:   PetscMPIInt    size;
146:   IS             istmp,nistmp;

150:   PetscLayoutGetLocalSize(is->map, &n);
151:   MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
152:   if (size == 1) {
153:     PetscMalloc1(n,&ii);
154:     for (i=0; i<n; i++) ii[idx[i]] = i;
155:     ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_OWN_POINTER,isout);
156:     ISSetPermutation(*isout);
157:   } else {
158:     /* crude, nonscalable get entire IS on each processor */
159:     ISAllGather(is,&istmp);
160:     ISSetPermutation(istmp);
161:     ISInvertPermutation(istmp,PETSC_DECIDE,&nistmp);
162:     ISDestroy(&istmp);
163:     /* get the part we need */
164:     if (nlocal == PETSC_DECIDE) nlocal = n;
165:     MPI_Scan(&nlocal,&nstart,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)is));
166:     if (PetscDefined(USE_DEBUG)) {
167:       PetscInt    N;
168:       PetscMPIInt rank;
169:       MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
170:       PetscLayoutGetSize(is->map, &N);
171:       if (rank == size-1) {
172:         if (nstart != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Sum of nlocal lengths %d != total IS length %d",nstart,N);
173:       }
174:     }
175:     nstart -= nlocal;
176:     ISGetIndices(nistmp,&idx);
177:     ISCreateGeneral(PetscObjectComm((PetscObject)is),nlocal,idx+nstart,PETSC_COPY_VALUES,isout);
178:     ISRestoreIndices(nistmp,&idx);
179:     ISDestroy(&nistmp);
180:   }
181:   return(0);
182: }

184: #if defined(PETSC_HAVE_HDF5)
185: static PetscErrorCode ISView_General_HDF5(IS is, PetscViewer viewer)
186: {
187:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
188:   hid_t           filespace;  /* file dataspace identifier */
189:   hid_t           chunkspace; /* chunk dataset property identifier */
190:   hid_t           dset_id;    /* dataset identifier */
191:   hid_t           memspace;   /* memory dataspace identifier */
192:   hid_t           inttype;    /* int type (H5T_NATIVE_INT or H5T_NATIVE_LLONG) */
193:   hid_t           file_id, group;
194:   hsize_t         dim, maxDims[3], dims[3], chunkDims[3], count[3],offset[3];
195:   PetscBool       timestepping;
196:   PetscInt        bs, N, n, timestep=PETSC_MIN_INT, low;
197:   hsize_t         chunksize;
198:   const PetscInt *ind;
199:   const char     *isname;
200:   PetscErrorCode  ierr;

203:   ISGetBlockSize(is,&bs);
204:   bs   = PetscMax(bs, 1); /* If N = 0, bs  = 0 as well */
205:   PetscViewerHDF5OpenGroup(viewer, &file_id, &group);
206:   PetscViewerHDF5IsTimestepping(viewer, &timestepping);
207:   if (timestepping) {
208:     PetscViewerHDF5GetTimestep(viewer, &timestep);
209:   }

211:   /* Create the dataspace for the dataset.
212:    *
213:    * dims - holds the current dimensions of the dataset
214:    *
215:    * maxDims - holds the maximum dimensions of the dataset (unlimited
216:    * for the number of time steps with the current dimensions for the
217:    * other dimensions; so only additional time steps can be added).
218:    *
219:    * chunkDims - holds the size of a single time step (required to
220:    * permit extending dataset).
221:    */
222:   dim = 0;
223:   chunksize = 1;
224:   if (timestep >= 0) {
225:     dims[dim]      = timestep+1;
226:     maxDims[dim]   = H5S_UNLIMITED;
227:     chunkDims[dim] = 1;
228:     ++dim;
229:   }
230:   ISGetSize(is, &N);
231:   ISGetLocalSize(is, &n);
232:   PetscHDF5IntCast(N/bs,dims + dim);

234:   maxDims[dim]   = dims[dim];
235:   chunkDims[dim] = PetscMax(1,dims[dim]);
236:   chunksize      *= chunkDims[dim];
237:   ++dim;
238:   if (bs >= 1) {
239:     dims[dim]      = bs;
240:     maxDims[dim]   = dims[dim];
241:     chunkDims[dim] = dims[dim];
242:     chunksize      *= chunkDims[dim];
243:     ++dim;
244:   }
245:   /* hdf5 chunks must be less than 4GB */
246:   if (chunksize > PETSC_HDF5_MAX_CHUNKSIZE/64) {
247:     if (bs >= 1) {
248:       if (chunkDims[dim-2] > (PetscInt)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE/64))) {
249:         chunkDims[dim-2] = (PetscInt)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE/64));
250:       } if (chunkDims[dim-1] > (PetscInt)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE/64))) {
251:         chunkDims[dim-1] = (PetscInt)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE/64));
252:       }
253:     } else {
254:       chunkDims[dim-1] = PETSC_HDF5_MAX_CHUNKSIZE/64;
255:     }
256:   }
257:   PetscStackCallHDF5Return(filespace,H5Screate_simple,(dim, dims, maxDims));

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

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

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

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

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

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

325:   /* Close/release resources */
326:   PetscStackCallHDF5(H5Gclose,(group));
327:   PetscStackCallHDF5(H5Sclose,(filespace));
328:   PetscStackCallHDF5(H5Sclose,(memspace));
329:   PetscStackCallHDF5(H5Dclose,(dset_id));

331:   PetscViewerHDF5WriteObjectAttribute(viewer,(PetscObject)is,"timestepping",PETSC_BOOL,&timestepping);
332:   PetscInfo1(is, "Wrote IS object with name %s\n", isname);
333:   return(0);
334: }
335: #endif

337: static PetscErrorCode ISView_General(IS is,PetscViewer viewer)
338: {
339:   IS_General     *sub = (IS_General*)is->data;
341:   PetscInt       i,n,*idx = sub->idx;
342:   PetscBool      iascii,isbinary,ishdf5;

345:   PetscLayoutGetLocalSize(is->map, &n);
346:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
347:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
348:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
349:   if (iascii) {
350:     MPI_Comm          comm;
351:     PetscMPIInt       rank,size;
352:     PetscViewerFormat fmt;
353:     PetscBool         isperm;

355:     PetscObjectGetComm((PetscObject)viewer,&comm);
356:     MPI_Comm_rank(comm,&rank);
357:     MPI_Comm_size(comm,&size);

359:     PetscViewerGetFormat(viewer,&fmt);
360:     ISGetInfo(is,IS_PERMUTATION,IS_LOCAL,PETSC_FALSE,&isperm);
361:     if (isperm && fmt != PETSC_VIEWER_ASCII_MATLAB) {PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");}
362:     PetscViewerASCIIPushSynchronized(viewer);
363:     if (size > 1) {
364:       if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
365:         const char* name;

367:         PetscObjectGetName((PetscObject)is,&name);
368:         PetscViewerASCIISynchronizedPrintf(viewer,"%s_%d = [...\n",name,rank);
369:         for (i=0; i<n; i++) {
370:           PetscViewerASCIISynchronizedPrintf(viewer,"%D\n",idx[i]+1);
371:         }
372:         PetscViewerASCIISynchronizedPrintf(viewer,"];\n");
373:       } else {
374:         PetscInt  st = 0;

376:         if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
377:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in set %D\n",rank,n);
378:         for (i=0; i<n; i++) {
379:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i + st,idx[i]);
380:         }
381:       }
382:     } else {
383:       if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
384:         const char* name;

386:         PetscObjectGetName((PetscObject)is,&name);
387:         PetscViewerASCIISynchronizedPrintf(viewer,"%s = [...\n",name);
388:         for (i=0; i<n; i++) {
389:           PetscViewerASCIISynchronizedPrintf(viewer,"%D\n",idx[i]+1);
390:         }
391:         PetscViewerASCIISynchronizedPrintf(viewer,"];\n");
392:       } else {
393:         PetscInt  st = 0;

395:         if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
396:         PetscViewerASCIISynchronizedPrintf(viewer,"Number of indices in set %D\n",n);
397:         for (i=0; i<n; i++) {
398:           PetscViewerASCIISynchronizedPrintf(viewer,"%D %D\n",i + st,idx[i]);
399:         }
400:       }
401:     }
402:     PetscViewerFlush(viewer);
403:     PetscViewerASCIIPopSynchronized(viewer);
404:   } else if (isbinary) {
405:     ISView_Binary(is,viewer);
406:   } else if (ishdf5) {
407: #if defined(PETSC_HAVE_HDF5)
408:     ISView_General_HDF5(is,viewer);
409: #endif
410:   }
411:   return(0);
412: }

414: static PetscErrorCode ISSort_General(IS is)
415: {
416:   IS_General     *sub = (IS_General*)is->data;
417:   PetscInt       n;

421:   PetscLayoutGetLocalSize(is->map, &n);
422:   PetscIntSortSemiOrdered(n,sub->idx);
423:   return(0);
424: }

426: static PetscErrorCode ISSortRemoveDups_General(IS is)
427: {
428:   IS_General     *sub = (IS_General*)is->data;
429:   PetscLayout    map;
430:   PetscInt       n;
431:   PetscBool      sorted;

435:   PetscLayoutGetLocalSize(is->map, &n);
436:   ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sorted);
437:   if (sorted) {
438:     PetscSortedRemoveDupsInt(&n,sub->idx);
439:   } else {
440:     PetscSortRemoveDupsInt(&n,sub->idx);
441:   }
442:   PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map);
443:   PetscLayoutDestroy(&is->map);
444:   is->map = map;
445:   return(0);
446: }

448: static PetscErrorCode ISSorted_General(IS is,PetscBool  *flg)
449: {

453:   ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,flg);
454:   return(0);
455: }

457: PetscErrorCode  ISToGeneral_General(IS is)
458: {
460:   return(0);
461: }

463: static struct _ISOps myops = { ISGetIndices_General,
464:                                ISRestoreIndices_General,
465:                                ISInvertPermutation_General,
466:                                ISSort_General,
467:                                ISSortRemoveDups_General,
468:                                ISSorted_General,
469:                                ISDuplicate_General,
470:                                ISDestroy_General,
471:                                ISView_General,
472:                                ISLoad_Default,
473:                                ISCopy_General,
474:                                ISToGeneral_General,
475:                                ISOnComm_General,
476:                                ISSetBlockSize_General,
477:                                ISContiguousLocal_General,
478:                                ISLocate_General,
479:                                /* no specializations of {sorted,unique,perm,interval}{local,global}
480:                                 * because the default checks in ISGetInfo_XXX in index.c are exactly
481:                                 * what we would do for ISGeneral */
482:                                NULL,
483:                                NULL,
484:                                NULL,
485:                                NULL,
486:                                NULL,
487:                                NULL,
488:                                NULL,
489:                                NULL};

491: PETSC_INTERN PetscErrorCode ISSetUp_General(IS);

493: PetscErrorCode ISSetUp_General(IS is)
494: {
496:   IS_General     *sub = (IS_General*)is->data;
497:   const PetscInt *idx = sub->idx;
498:   PetscInt       n,i,min,max;

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

503:   if (n) {
504:     min = max = idx[0];
505:     for (i=1; i<n; i++) {
506:       if (idx[i] < min) min = idx[i];
507:       if (idx[i] > max) max = idx[i];
508:     }
509:     is->min = min;
510:     is->max = max;
511:   } else {
512:     is->min = PETSC_MAX_INT;
513:     is->max = PETSC_MIN_INT;
514:   }
515:   return(0);
516: }

518: /*@
519:    ISCreateGeneral - Creates a data structure for an index set
520:    containing a list of integers.

522:    Collective

524:    Input Parameters:
525: +  comm - the MPI communicator
526: .  n - the length of the index set
527: .  idx - the list of integers
528: -  mode - PETSC_COPY_VALUES, PETSC_OWN_POINTER, or PETSC_USE_POINTER; see PetscCopyMode for meaning of this flag.

530:    Output Parameter:
531: .  is - the new index set

533:    Notes:
534:    When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
535:    conceptually the same as MPI_Group operations. The IS are then
536:    distributed sets of indices and thus certain operations on them are
537:    collective.

539:    Level: beginner

541: .seealso: ISCreateStride(), ISCreateBlock(), ISAllGather(), PETSC_COPY_VALUES, PETSC_OWN_POINTER, PETSC_USE_POINTER, PetscCopyMode
542: @*/
543: PetscErrorCode  ISCreateGeneral(MPI_Comm comm,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
544: {

548:   ISCreate(comm,is);
549:   ISSetType(*is,ISGENERAL);
550:   ISGeneralSetIndices(*is,n,idx,mode);
551:   return(0);
552: }

554: /*@
555:    ISGeneralSetIndices - Sets the indices for an ISGENERAL index set

557:    Collective on IS

559:    Input Parameters:
560: +  is - the index set
561: .  n - the length of the index set
562: .  idx - the list of integers
563: -  mode - see PetscCopyMode for meaning of this flag.

565:    Level: beginner

567: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather(), ISBlockSetIndices(), ISGENERAL, PetscCopyMode
568: @*/
569: PetscErrorCode  ISGeneralSetIndices(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
570: {

574:   ISClearInfoCache(is,PETSC_FALSE);
575:   PetscUseMethod(is,"ISGeneralSetIndices_C",(IS,PetscInt,const PetscInt[],PetscCopyMode),(is,n,idx,mode));
576:   return(0);
577: }

579: PetscErrorCode  ISGeneralSetIndices_General(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
580: {
581:   PetscLayout    map;
583:   IS_General     *sub = (IS_General*)is->data;

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

589:   PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is),n,PETSC_DECIDE,is->map->bs,&map);
590:   PetscLayoutDestroy(&is->map);
591:   is->map = map;

593:   if (sub->allocated) {PetscFree(sub->idx);}
594:   if (mode == PETSC_COPY_VALUES) {
595:     PetscMalloc1(n,&sub->idx);
596:     PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
597:     PetscArraycpy(sub->idx,idx,n);
598:     sub->allocated = PETSC_TRUE;
599:   } else if (mode == PETSC_OWN_POINTER) {
600:     sub->idx = (PetscInt*)idx;
601:     PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
602:     sub->allocated = PETSC_TRUE;
603:   } else {
604:     sub->idx = (PetscInt*)idx;
605:     sub->allocated = PETSC_FALSE;
606:   }

608:   ISSetUp_General(is);
609:   ISViewFromOptions(is,NULL,"-is_view");
610:   return(0);
611: }

613: static PetscErrorCode ISGeneralFilter_General(IS is, PetscInt start, PetscInt end)
614: {
615:   IS_General     *sub = (IS_General*)is->data;
616:   PetscInt       *idx = sub->idx,*idxnew;
617:   PetscInt       i,n = is->map->n,nnew = 0,o;

621:   for (i=0; i<n; ++i)
622:     if (idx[i] >= start && idx[i] < end)
623:       nnew++;
624:   PetscMalloc1(nnew, &idxnew);
625:   for (o=0, i=0; i<n; i++) {
626:     if (idx[i] >= start && idx[i] < end)
627:       idxnew[o++] = idx[i];
628:   }
629:   ISGeneralSetIndices_General(is,nnew,idxnew,PETSC_OWN_POINTER);
630:   return(0);
631: }

633: /*@
634:    ISGeneralFilter - Remove all indices outside of [start, end)

636:    Collective on IS

638:    Input Parameters:
639: +  is - the index set
640: .  start - the lowest index kept
641: -  end - one more than the highest index kept

643:    Level: beginner

645: .seealso: ISCreateGeneral(), ISGeneralSetIndices()
646: @*/
647: PetscErrorCode ISGeneralFilter(IS is, PetscInt start, PetscInt end)
648: {

653:   ISClearInfoCache(is,PETSC_FALSE);
654:   PetscUseMethod(is,"ISGeneralFilter_C",(IS,PetscInt,PetscInt),(is,start,end));
655:   return(0);
656: }

658: PETSC_EXTERN PetscErrorCode ISCreate_General(IS is)
659: {
661:   IS_General     *sub;

664:   PetscNewLog(is,&sub);
665:   is->data = (void *) sub;
666:   PetscMemcpy(is->ops,&myops,sizeof(myops));
667:   PetscObjectComposeFunction((PetscObject)is,"ISGeneralSetIndices_C",ISGeneralSetIndices_General);
668:   PetscObjectComposeFunction((PetscObject)is,"ISGeneralFilter_C",ISGeneralFilter_General);
669:   return(0);
670: }