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: {
10: IS_General *sub = (IS_General*)is->data;
11: PetscInt n;
13: PetscLayoutGetLocalSize(is->map, &n);
14: ISCreateGeneral(PetscObjectComm((PetscObject) is), n, sub->idx, PETSC_COPY_VALUES, newIS);
15: return 0;
16: }
18: static PetscErrorCode ISDestroy_General(IS is)
19: {
20: IS_General *is_general = (IS_General*)is->data;
22: if (is_general->allocated) PetscFree(is_general->idx);
23: PetscObjectComposeFunction((PetscObject)is,"ISGeneralSetIndices_C",NULL);
24: PetscObjectComposeFunction((PetscObject)is,"ISGeneralFilter_C",NULL);
25: PetscFree(is->data);
26: return 0;
27: }
29: static PetscErrorCode ISCopy_General(IS is,IS isy)
30: {
31: IS_General *is_general = (IS_General*)is->data,*isy_general = (IS_General*)isy->data;
32: PetscInt n, N, ny, Ny;
34: PetscLayoutGetLocalSize(is->map, &n);
35: PetscLayoutGetSize(is->map, &N);
36: PetscLayoutGetLocalSize(isy->map, &ny);
37: PetscLayoutGetSize(isy->map, &Ny);
39: PetscArraycpy(isy_general->idx,is_general->idx,n);
40: return 0;
41: }
43: static PetscErrorCode ISOnComm_General(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
44: {
45: IS_General *sub = (IS_General*)is->data;
46: PetscInt n;
49: PetscLayoutGetLocalSize(is->map, &n);
50: ISCreateGeneral(comm,n,sub->idx,mode,newis);
51: return 0;
52: }
54: static PetscErrorCode ISSetBlockSize_General(IS is,PetscInt bs)
55: {
56: PetscLayoutSetBlockSize(is->map, bs);
57: return 0;
58: }
60: static PetscErrorCode ISContiguousLocal_General(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
61: {
62: IS_General *sub = (IS_General*)is->data;
63: PetscInt n,i,p;
65: *start = 0;
66: *contig = PETSC_TRUE;
67: PetscLayoutGetLocalSize(is->map, &n);
68: if (!n) return 0;
69: p = sub->idx[0];
70: if (p < gstart) goto nomatch;
71: *start = p - gstart;
72: if (n > gend-p) goto nomatch;
73: for (i=1; i<n; i++,p++) {
74: if (sub->idx[i] != p+1) goto nomatch;
75: }
76: return 0;
77: nomatch:
78: *start = -1;
79: *contig = PETSC_FALSE;
80: return 0;
81: }
83: static PetscErrorCode ISLocate_General(IS is,PetscInt key,PetscInt *location)
84: {
85: IS_General *sub = (IS_General*)is->data;
86: PetscInt numIdx, i;
87: PetscBool sorted;
89: PetscLayoutGetLocalSize(is->map,&numIdx);
90: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sorted);
91: if (sorted) PetscFindInt(key,numIdx,sub->idx,location);
92: else {
93: const PetscInt *idx = sub->idx;
95: *location = -1;
96: for (i = 0; i < numIdx; i++) {
97: if (idx[i] == key) {
98: *location = i;
99: return 0;
100: }
101: }
102: }
103: return 0;
104: }
106: static PetscErrorCode ISGetIndices_General(IS in,const PetscInt *idx[])
107: {
108: IS_General *sub = (IS_General*)in->data;
110: *idx = sub->idx;
111: return 0;
112: }
114: static PetscErrorCode ISRestoreIndices_General(IS in,const PetscInt *idx[])
115: {
116: IS_General *sub = (IS_General*)in->data;
118: /* F90Array1dCreate() inside ISRestoreArrayF90() does not keep array when zero length array */
120: return 0;
121: }
123: static PetscErrorCode ISInvertPermutation_General(IS is,PetscInt nlocal,IS *isout)
124: {
125: IS_General *sub = (IS_General*)is->data;
126: PetscInt i,*ii,n,nstart;
127: const PetscInt *idx = sub->idx;
128: PetscMPIInt size;
129: IS istmp,nistmp;
131: PetscLayoutGetLocalSize(is->map, &n);
132: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
133: if (size == 1) {
134: PetscMalloc1(n,&ii);
135: for (i=0; i<n; i++) ii[idx[i]] = i;
136: ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_OWN_POINTER,isout);
137: ISSetPermutation(*isout);
138: } else {
139: /* crude, nonscalable get entire IS on each processor */
140: ISAllGather(is,&istmp);
141: ISSetPermutation(istmp);
142: ISInvertPermutation(istmp,PETSC_DECIDE,&nistmp);
143: ISDestroy(&istmp);
144: /* get the part we need */
145: if (nlocal == PETSC_DECIDE) nlocal = n;
146: MPI_Scan(&nlocal,&nstart,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)is));
147: if (PetscDefined(USE_DEBUG)) {
148: PetscInt N;
149: PetscMPIInt rank;
150: MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
151: PetscLayoutGetSize(is->map, &N);
153: }
154: nstart -= nlocal;
155: ISGetIndices(nistmp,&idx);
156: ISCreateGeneral(PetscObjectComm((PetscObject)is),nlocal,idx+nstart,PETSC_COPY_VALUES,isout);
157: ISRestoreIndices(nistmp,&idx);
158: ISDestroy(&nistmp);
159: }
160: return 0;
161: }
163: #if defined(PETSC_HAVE_HDF5)
164: static PetscErrorCode ISView_General_HDF5(IS is, PetscViewer viewer)
165: {
166: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
167: hid_t filespace; /* file dataspace identifier */
168: hid_t chunkspace; /* chunk dataset property identifier */
169: hid_t dset_id; /* dataset identifier */
170: hid_t memspace; /* memory dataspace identifier */
171: hid_t inttype; /* int type (H5T_NATIVE_INT or H5T_NATIVE_LLONG) */
172: hid_t file_id, group;
173: hsize_t dim, maxDims[3], dims[3], chunkDims[3], count[3],offset[3];
174: PetscBool timestepping;
175: PetscInt bs, N, n, timestep=PETSC_MIN_INT, low;
176: hsize_t chunksize;
177: const PetscInt *ind;
178: const char *isname;
180: ISGetBlockSize(is,&bs);
181: bs = PetscMax(bs, 1); /* If N = 0, bs = 0 as well */
182: PetscViewerHDF5OpenGroup(viewer, &file_id, &group);
183: PetscViewerHDF5IsTimestepping(viewer, ×tepping);
184: if (timestepping) {
185: PetscViewerHDF5GetTimestep(viewer, ×tep);
186: }
188: /* Create the dataspace for the dataset.
189: *
190: * dims - holds the current dimensions of the dataset
191: *
192: * maxDims - holds the maximum dimensions of the dataset (unlimited
193: * for the number of time steps with the current dimensions for the
194: * other dimensions; so only additional time steps can be added).
195: *
196: * chunkDims - holds the size of a single time step (required to
197: * permit extending dataset).
198: */
199: dim = 0;
200: chunksize = 1;
201: if (timestep >= 0) {
202: dims[dim] = timestep+1;
203: maxDims[dim] = H5S_UNLIMITED;
204: chunkDims[dim] = 1;
205: ++dim;
206: }
207: ISGetSize(is, &N);
208: ISGetLocalSize(is, &n);
209: PetscHDF5IntCast(N/bs,dims + dim);
211: maxDims[dim] = dims[dim];
212: chunkDims[dim] = PetscMax(1,dims[dim]);
213: chunksize *= chunkDims[dim];
214: ++dim;
215: if (bs >= 1) {
216: dims[dim] = bs;
217: maxDims[dim] = dims[dim];
218: chunkDims[dim] = dims[dim];
219: chunksize *= chunkDims[dim];
220: ++dim;
221: }
222: /* hdf5 chunks must be less than 4GB */
223: if (chunksize > PETSC_HDF5_MAX_CHUNKSIZE/64) {
224: if (bs >= 1) {
225: if (chunkDims[dim-2] > (PetscInt)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE/64))) {
226: chunkDims[dim-2] = (PetscInt)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE/64));
227: } if (chunkDims[dim-1] > (PetscInt)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE/64))) {
228: chunkDims[dim-1] = (PetscInt)PetscSqrtReal((PetscReal)(PETSC_HDF5_MAX_CHUNKSIZE/64));
229: }
230: } else {
231: chunkDims[dim-1] = PETSC_HDF5_MAX_CHUNKSIZE/64;
232: }
233: }
234: PetscStackCallHDF5Return(filespace,H5Screate_simple,(dim, dims, maxDims));
236: #if defined(PETSC_USE_64BIT_INDICES)
237: inttype = H5T_NATIVE_LLONG;
238: #else
239: inttype = H5T_NATIVE_INT;
240: #endif
242: /* Create the dataset with default properties and close filespace */
243: PetscObjectGetName((PetscObject) is, &isname);
244: if (!H5Lexists(group, isname, H5P_DEFAULT)) {
245: /* Create chunk */
246: PetscStackCallHDF5Return(chunkspace,H5Pcreate,(H5P_DATASET_CREATE));
247: PetscStackCallHDF5(H5Pset_chunk,(chunkspace, dim, chunkDims));
249: PetscStackCallHDF5Return(dset_id,H5Dcreate2,(group, isname, inttype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
250: PetscStackCallHDF5(H5Pclose,(chunkspace));
251: } else {
252: PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, isname, H5P_DEFAULT));
253: PetscStackCallHDF5(H5Dset_extent,(dset_id, dims));
254: }
255: PetscStackCallHDF5(H5Sclose,(filespace));
257: /* Each process defines a dataset and writes it to the hyperslab in the file */
258: dim = 0;
259: if (timestep >= 0) {
260: count[dim] = 1;
261: ++dim;
262: }
263: PetscHDF5IntCast(n/bs,count + dim);
264: ++dim;
265: if (bs >= 1) {
266: count[dim] = bs;
267: ++dim;
268: }
269: if (n > 0 || H5_VERSION_GE(1,10,0)) {
270: PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL));
271: } else {
272: /* Can't create dataspace with zero for any dimension, so create null dataspace. */
273: PetscStackCallHDF5Return(memspace,H5Screate,(H5S_NULL));
274: }
276: /* Select hyperslab in the file */
277: PetscLayoutGetRange(is->map, &low, NULL);
278: dim = 0;
279: if (timestep >= 0) {
280: offset[dim] = timestep;
281: ++dim;
282: }
283: PetscHDF5IntCast(low/bs,offset + dim);
284: ++dim;
285: if (bs >= 1) {
286: offset[dim] = 0;
287: ++dim;
288: }
289: if (n > 0 || H5_VERSION_GE(1,10,0)) {
290: PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id));
291: PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
292: } else {
293: /* Create null filespace to match null memspace. */
294: PetscStackCallHDF5Return(filespace,H5Screate,(H5S_NULL));
295: }
297: ISGetIndices(is, &ind);
298: PetscStackCallHDF5(H5Dwrite,(dset_id, inttype, memspace, filespace, hdf5->dxpl_id, ind));
299: PetscStackCallHDF5(H5Fflush,(file_id, H5F_SCOPE_GLOBAL));
300: ISRestoreIndices(is, &ind);
302: /* Close/release resources */
303: PetscStackCallHDF5(H5Gclose,(group));
304: PetscStackCallHDF5(H5Sclose,(filespace));
305: PetscStackCallHDF5(H5Sclose,(memspace));
306: PetscStackCallHDF5(H5Dclose,(dset_id));
308: if (timestepping) {
309: PetscViewerHDF5WriteObjectAttribute(viewer,(PetscObject)is,"timestepping",PETSC_BOOL,×tepping);
310: }
311: PetscInfo(is, "Wrote IS object with name %s\n", isname);
312: return 0;
313: }
314: #endif
316: static PetscErrorCode ISView_General(IS is,PetscViewer viewer)
317: {
318: IS_General *sub = (IS_General*)is->data;
319: PetscInt i,n,*idx = sub->idx;
320: PetscBool iascii,isbinary,ishdf5;
322: PetscLayoutGetLocalSize(is->map, &n);
323: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
324: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
325: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
326: if (iascii) {
327: MPI_Comm comm;
328: PetscMPIInt rank,size;
329: PetscViewerFormat fmt;
330: PetscBool isperm;
332: PetscObjectGetComm((PetscObject)viewer,&comm);
333: MPI_Comm_rank(comm,&rank);
334: MPI_Comm_size(comm,&size);
336: PetscViewerGetFormat(viewer,&fmt);
337: ISGetInfo(is,IS_PERMUTATION,IS_LOCAL,PETSC_FALSE,&isperm);
338: if (isperm && fmt != PETSC_VIEWER_ASCII_MATLAB) PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");
339: PetscViewerASCIIPushSynchronized(viewer);
340: if (size > 1) {
341: if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
342: const char* name;
344: PetscObjectGetName((PetscObject)is,&name);
345: PetscViewerASCIISynchronizedPrintf(viewer,"%s_%d = [...\n",name,rank);
346: for (i=0; i<n; i++) {
347: PetscViewerASCIISynchronizedPrintf(viewer,"%" PetscInt_FMT "\n",idx[i]+1);
348: }
349: PetscViewerASCIISynchronizedPrintf(viewer,"];\n");
350: } else {
351: PetscInt st = 0;
353: if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
354: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in set %" PetscInt_FMT "\n",rank,n);
355: for (i=0; i<n; i++) {
356: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %" PetscInt_FMT " %" PetscInt_FMT "\n",rank,i + st,idx[i]);
357: }
358: }
359: } else {
360: if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
361: const char* name;
363: PetscObjectGetName((PetscObject)is,&name);
364: PetscViewerASCIISynchronizedPrintf(viewer,"%s = [...\n",name);
365: for (i=0; i<n; i++) {
366: PetscViewerASCIISynchronizedPrintf(viewer,"%" PetscInt_FMT "\n",idx[i]+1);
367: }
368: PetscViewerASCIISynchronizedPrintf(viewer,"];\n");
369: } else {
370: PetscInt st = 0;
372: if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
373: PetscViewerASCIISynchronizedPrintf(viewer,"Number of indices in set %" PetscInt_FMT "\n",n);
374: for (i=0; i<n; i++) {
375: PetscViewerASCIISynchronizedPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT "\n",i + st,idx[i]);
376: }
377: }
378: }
379: PetscViewerFlush(viewer);
380: PetscViewerASCIIPopSynchronized(viewer);
381: } else if (isbinary) {
382: ISView_Binary(is,viewer);
383: } else if (ishdf5) {
384: #if defined(PETSC_HAVE_HDF5)
385: ISView_General_HDF5(is,viewer);
386: #endif
387: }
388: return 0;
389: }
391: static PetscErrorCode ISSort_General(IS is)
392: {
393: IS_General *sub = (IS_General*)is->data;
394: PetscInt n;
396: PetscLayoutGetLocalSize(is->map, &n);
397: PetscIntSortSemiOrdered(n,sub->idx);
398: return 0;
399: }
401: static PetscErrorCode ISSortRemoveDups_General(IS is)
402: {
403: IS_General *sub = (IS_General*)is->data;
404: PetscLayout map;
405: PetscInt n;
406: PetscBool sorted;
408: PetscLayoutGetLocalSize(is->map, &n);
409: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sorted);
410: if (sorted) {
411: PetscSortedRemoveDupsInt(&n,sub->idx);
412: } else {
413: PetscSortRemoveDupsInt(&n,sub->idx);
414: }
415: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map);
416: PetscLayoutDestroy(&is->map);
417: is->map = map;
418: return 0;
419: }
421: static PetscErrorCode ISSorted_General(IS is,PetscBool *flg)
422: {
423: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,flg);
424: return 0;
425: }
427: PetscErrorCode ISToGeneral_General(IS is)
428: {
429: return 0;
430: }
432: static struct _ISOps myops = { ISGetIndices_General,
433: ISRestoreIndices_General,
434: ISInvertPermutation_General,
435: ISSort_General,
436: ISSortRemoveDups_General,
437: ISSorted_General,
438: ISDuplicate_General,
439: ISDestroy_General,
440: ISView_General,
441: ISLoad_Default,
442: ISCopy_General,
443: ISToGeneral_General,
444: ISOnComm_General,
445: ISSetBlockSize_General,
446: ISContiguousLocal_General,
447: ISLocate_General,
448: /* no specializations of {sorted,unique,perm,interval}{local,global}
449: * because the default checks in ISGetInfo_XXX in index.c are exactly
450: * what we would do for ISGeneral */
451: NULL,
452: NULL,
453: NULL,
454: NULL,
455: NULL,
456: NULL,
457: NULL,
458: NULL};
460: PETSC_INTERN PetscErrorCode ISSetUp_General(IS);
462: PetscErrorCode ISSetUp_General(IS is)
463: {
464: IS_General *sub = (IS_General*)is->data;
465: const PetscInt *idx = sub->idx;
466: PetscInt n,i,min,max;
468: PetscLayoutGetLocalSize(is->map, &n);
470: if (n) {
471: min = max = idx[0];
472: for (i=1; i<n; i++) {
473: if (idx[i] < min) min = idx[i];
474: if (idx[i] > max) max = idx[i];
475: }
476: is->min = min;
477: is->max = max;
478: } else {
479: is->min = PETSC_MAX_INT;
480: is->max = PETSC_MIN_INT;
481: }
482: return 0;
483: }
485: /*@
486: ISCreateGeneral - Creates a data structure for an index set
487: containing a list of integers.
489: Collective
491: Input Parameters:
492: + comm - the MPI communicator
493: . n - the length of the index set
494: . idx - the list of integers
495: - mode - PETSC_COPY_VALUES, PETSC_OWN_POINTER, or PETSC_USE_POINTER; see PetscCopyMode for meaning of this flag.
497: Output Parameter:
498: . is - the new index set
500: Notes:
501: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
502: conceptually the same as MPI_Group operations. The IS are then
503: distributed sets of indices and thus certain operations on them are
504: collective.
506: Level: beginner
508: .seealso: ISCreateStride(), ISCreateBlock(), ISAllGather(), PETSC_COPY_VALUES, PETSC_OWN_POINTER, PETSC_USE_POINTER, PetscCopyMode
509: @*/
510: PetscErrorCode ISCreateGeneral(MPI_Comm comm,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
511: {
512: ISCreate(comm,is);
513: ISSetType(*is,ISGENERAL);
514: ISGeneralSetIndices(*is,n,idx,mode);
515: return 0;
516: }
518: /*@
519: ISGeneralSetIndices - Sets the indices for an ISGENERAL index set
521: Collective on IS
523: Input Parameters:
524: + is - the index set
525: . n - the length of the index set
526: . idx - the list of integers
527: - mode - see PetscCopyMode for meaning of this flag.
529: Level: beginner
531: .seealso: ISCreateGeneral(), ISGeneralSetIndicesFromMask(), ISBlockSetIndices(), ISGENERAL, PetscCopyMode
532: @*/
533: PetscErrorCode ISGeneralSetIndices(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
534: {
537: ISClearInfoCache(is,PETSC_FALSE);
538: PetscUseMethod(is,"ISGeneralSetIndices_C",(IS,PetscInt,const PetscInt[],PetscCopyMode),(is,n,idx,mode));
539: return 0;
540: }
542: PetscErrorCode ISGeneralSetIndices_General(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
543: {
544: PetscLayout map;
545: IS_General *sub = (IS_General*)is->data;
550: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is),n,PETSC_DECIDE,is->map->bs,&map);
551: PetscLayoutDestroy(&is->map);
552: is->map = map;
554: if (sub->allocated) PetscFree(sub->idx);
555: if (mode == PETSC_COPY_VALUES) {
556: PetscMalloc1(n,&sub->idx);
557: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
558: PetscArraycpy(sub->idx,idx,n);
559: sub->allocated = PETSC_TRUE;
560: } else if (mode == PETSC_OWN_POINTER) {
561: sub->idx = (PetscInt*)idx;
562: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
563: sub->allocated = PETSC_TRUE;
564: } else {
565: sub->idx = (PetscInt*)idx;
566: sub->allocated = PETSC_FALSE;
567: }
569: ISSetUp_General(is);
570: ISViewFromOptions(is,NULL,"-is_view");
571: return 0;
572: }
574: /*@
575: ISGeneralSetIndicesFromMask - Sets the indices for an ISGENERAL index set using a boolean mask
577: Collective on IS
579: Input Parameters:
580: + is - the index set
581: . rstart - the range start index (inclusive)
582: . rend - the range end index (exclusive)
583: - mask - the boolean mask array of length rend-rstart, indices will be set for each PETSC_TRUE value in the array
585: Notes:
586: The mask array may be freed by the user after this call.
588: Example:
589: $ PetscBool mask[] = {PETSC_FALSE, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE, PETSC_TRUE};
590: $ ISGeneralSetIndicesFromMask(is,10,15,mask);
591: will feed the IS with indices
592: $ {11, 14}
593: locally.
595: Level: beginner
597: .seealso: ISCreateGeneral(), ISGeneralSetIndices(), ISGENERAL
598: @*/
599: PetscErrorCode ISGeneralSetIndicesFromMask(IS is,PetscInt rstart,PetscInt rend,const PetscBool mask[])
600: {
603: ISClearInfoCache(is,PETSC_FALSE);
604: PetscUseMethod(is,"ISGeneralSetIndicesFromMask_C",(IS,PetscInt,PetscInt,const PetscBool[]),(is,rstart,rend,mask));
605: return 0;
606: }
608: PetscErrorCode ISGeneralSetIndicesFromMask_General(IS is,PetscInt rstart,PetscInt rend,const PetscBool mask[])
609: {
610: PetscInt i,nidx;
611: PetscInt *idx;
613: for (i=0,nidx=0; i<rend-rstart; i++) if (mask[i]) nidx++;
614: PetscMalloc1(nidx,&idx);
615: for (i=0,nidx=0; i<rend-rstart; i++) {
616: if (mask[i]) {
617: idx[nidx] = i+rstart;
618: nidx++;
619: }
620: }
621: ISGeneralSetIndices_General(is,nidx,idx,PETSC_OWN_POINTER);
622: return 0;
623: }
625: static PetscErrorCode ISGeneralFilter_General(IS is, PetscInt start, PetscInt end)
626: {
627: IS_General *sub = (IS_General*)is->data;
628: PetscInt *idx = sub->idx,*idxnew;
629: PetscInt i,n = is->map->n,nnew = 0,o;
631: for (i=0; i<n; ++i)
632: if (idx[i] >= start && idx[i] < end)
633: nnew++;
634: PetscMalloc1(nnew, &idxnew);
635: for (o=0, i=0; i<n; i++) {
636: if (idx[i] >= start && idx[i] < end)
637: idxnew[o++] = idx[i];
638: }
639: ISGeneralSetIndices_General(is,nnew,idxnew,PETSC_OWN_POINTER);
640: return 0;
641: }
643: /*@
644: ISGeneralFilter - Remove all indices outside of [start, end)
646: Collective on IS
648: Input Parameters:
649: + is - the index set
650: . start - the lowest index kept
651: - end - one more than the highest index kept
653: Level: beginner
655: .seealso: ISCreateGeneral(), ISGeneralSetIndices()
656: @*/
657: PetscErrorCode ISGeneralFilter(IS is, PetscInt start, PetscInt end)
658: {
660: ISClearInfoCache(is,PETSC_FALSE);
661: PetscUseMethod(is,"ISGeneralFilter_C",(IS,PetscInt,PetscInt),(is,start,end));
662: return 0;
663: }
665: PETSC_EXTERN PetscErrorCode ISCreate_General(IS is)
666: {
667: IS_General *sub;
669: PetscNewLog(is,&sub);
670: is->data = (void *) sub;
671: PetscMemcpy(is->ops,&myops,sizeof(myops));
672: PetscObjectComposeFunction((PetscObject)is,"ISGeneralSetIndices_C",ISGeneralSetIndices_General);
673: PetscObjectComposeFunction((PetscObject)is,"ISGeneralSetIndicesFromMask_C",ISGeneralSetIndicesFromMask_General);
674: PetscObjectComposeFunction((PetscObject)is,"ISGeneralFilter_C",ISGeneralFilter_General);
675: return 0;
676: }