Actual source code: general.c
petsc-master 2020-02-15
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 ISCopy_General(IS is,IS isy)
35: {
36: IS_General *is_general = (IS_General*)is->data,*isy_general = (IS_General*)isy->data;
37: PetscInt n, N, ny, Ny;
41: PetscLayoutGetLocalSize(is->map, &n);
42: PetscLayoutGetSize(is->map, &N);
43: PetscLayoutGetLocalSize(isy->map, &ny);
44: PetscLayoutGetSize(isy->map, &Ny);
45: if (n != ny || N != Ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Index sets incompatible");
46: PetscArraycpy(isy_general->idx,is_general->idx,n);
47: return(0);
48: }
50: static PetscErrorCode ISOnComm_General(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
51: {
53: IS_General *sub = (IS_General*)is->data;
54: PetscInt n;
57: if (mode == PETSC_OWN_POINTER) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Cannot use PETSC_OWN_POINTER");
58: PetscLayoutGetLocalSize(is->map, &n);
59: ISCreateGeneral(comm,n,sub->idx,mode,newis);
60: return(0);
61: }
63: static PetscErrorCode ISSetBlockSize_General(IS is,PetscInt bs)
64: {
68: PetscLayoutSetBlockSize(is->map, bs);
69: return(0);
70: }
72: static PetscErrorCode ISContiguousLocal_General(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
73: {
74: IS_General *sub = (IS_General*)is->data;
75: PetscInt n,i,p;
79: *start = 0;
80: *contig = PETSC_TRUE;
81: PetscLayoutGetLocalSize(is->map, &n);
82: if (!n) return(0);
83: p = sub->idx[0];
84: if (p < gstart) goto nomatch;
85: *start = p - gstart;
86: if (n > gend-p) goto nomatch;
87: for (i=1; i<n; i++,p++) {
88: if (sub->idx[i] != p+1) goto nomatch;
89: }
90: return(0);
91: nomatch:
92: *start = -1;
93: *contig = PETSC_FALSE;
94: return(0);
95: }
97: static PetscErrorCode ISLocate_General(IS is,PetscInt key,PetscInt *location)
98: {
99: IS_General *sub = (IS_General*)is->data;
100: PetscInt numIdx, i;
101: PetscBool sorted;
105: PetscLayoutGetLocalSize(is->map,&numIdx);
106: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sorted);
107: if (sorted) {PetscFindInt(key,numIdx,sub->idx,location);}
108: else {
109: const PetscInt *idx = sub->idx;
111: *location = -1;
112: for (i = 0; i < numIdx; i++) {
113: if (idx[i] == key) {
114: *location = i;
115: return(0);
116: }
117: }
118: }
119: return(0);
120: }
122: static PetscErrorCode ISGetIndices_General(IS in,const PetscInt *idx[])
123: {
124: IS_General *sub = (IS_General*)in->data;
127: *idx = sub->idx;
128: return(0);
129: }
131: static PetscErrorCode ISRestoreIndices_General(IS in,const PetscInt *idx[])
132: {
133: IS_General *sub = (IS_General*)in->data;
136: if (*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 defined(PETSC_USE_DEBUG)
167: {
168: PetscInt N;
169: PetscMPIInt rank;
170: MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
171: PetscLayoutGetSize(is->map, &N);
172: if (rank == size-1) {
173: if (nstart != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Sum of nlocal lengths %d != total IS length %d",nstart,N);
174: }
175: }
176: #endif
177: nstart -= nlocal;
178: ISGetIndices(nistmp,&idx);
179: ISCreateGeneral(PetscObjectComm((PetscObject)is),nlocal,idx+nstart,PETSC_COPY_VALUES,isout);
180: ISRestoreIndices(nistmp,&idx);
181: ISDestroy(&nistmp);
182: }
183: return(0);
184: }
186: #if defined(PETSC_HAVE_HDF5)
187: static PetscErrorCode ISView_General_HDF5(IS is, PetscViewer viewer)
188: {
189: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
190: hid_t filespace; /* file dataspace identifier */
191: hid_t chunkspace; /* chunk dataset property identifier */
192: hid_t dset_id; /* dataset identifier */
193: hid_t memspace; /* memory dataspace identifier */
194: hid_t inttype; /* int type (H5T_NATIVE_INT or H5T_NATIVE_LLONG) */
195: hid_t file_id, group;
196: hsize_t dim, maxDims[3], dims[3], chunkDims[3], count[3],offset[3];
197: PetscInt bs, N, n, timestep, low;
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: PetscViewerHDF5GetTimestep(viewer, ×tep);
208: /* Create the dataspace for the dataset.
209: *
210: * dims - holds the current dimensions of the dataset
211: *
212: * maxDims - holds the maximum dimensions of the dataset (unlimited
213: * for the number of time steps with the current dimensions for the
214: * other dimensions; so only additional time steps can be added).
215: *
216: * chunkDims - holds the size of a single time step (required to
217: * permit extending dataset).
218: */
219: dim = 0;
220: if (timestep >= 0) {
221: dims[dim] = timestep+1;
222: maxDims[dim] = H5S_UNLIMITED;
223: chunkDims[dim] = 1;
224: ++dim;
225: }
226: ISGetSize(is, &N);
227: ISGetLocalSize(is, &n);
228: PetscHDF5IntCast(N/bs,dims + dim);
230: maxDims[dim] = dims[dim];
231: chunkDims[dim] = PetscMax(1,dims[dim]);
232: ++dim;
233: if (bs >= 1) {
234: dims[dim] = bs;
235: maxDims[dim] = dims[dim];
236: chunkDims[dim] = dims[dim];
237: ++dim;
238: }
239: PetscStackCallHDF5Return(filespace,H5Screate_simple,(dim, dims, maxDims));
241: #if defined(PETSC_USE_64BIT_INDICES)
242: inttype = H5T_NATIVE_LLONG;
243: #else
244: inttype = H5T_NATIVE_INT;
245: #endif
247: /* Create the dataset with default properties and close filespace */
248: PetscObjectGetName((PetscObject) is, &isname);
249: if (!H5Lexists(group, isname, H5P_DEFAULT)) {
250: /* Create chunk */
251: PetscStackCallHDF5Return(chunkspace,H5Pcreate,(H5P_DATASET_CREATE));
252: PetscStackCallHDF5(H5Pset_chunk,(chunkspace, dim, chunkDims));
254: PetscStackCallHDF5Return(dset_id,H5Dcreate2,(group, isname, inttype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
255: PetscStackCallHDF5(H5Pclose,(chunkspace));
256: } else {
257: PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, isname, H5P_DEFAULT));
258: PetscStackCallHDF5(H5Dset_extent,(dset_id, dims));
259: }
260: PetscStackCallHDF5(H5Sclose,(filespace));
262: /* Each process defines a dataset and writes it to the hyperslab in the file */
263: dim = 0;
264: if (timestep >= 0) {
265: count[dim] = 1;
266: ++dim;
267: }
268: PetscHDF5IntCast(n/bs,count + dim);
269: ++dim;
270: if (bs >= 1) {
271: count[dim] = bs;
272: ++dim;
273: }
274: if (n > 0) {
275: PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL));
276: } else {
277: /* Can't create dataspace with zero for any dimension, so create null dataspace. */
278: PetscStackCallHDF5Return(memspace,H5Screate,(H5S_NULL));
279: }
281: /* Select hyperslab in the file */
282: PetscLayoutGetRange(is->map, &low, NULL);
283: dim = 0;
284: if (timestep >= 0) {
285: offset[dim] = timestep;
286: ++dim;
287: }
288: PetscHDF5IntCast(low/bs,offset + dim);
289: ++dim;
290: if (bs >= 1) {
291: offset[dim] = 0;
292: ++dim;
293: }
294: if (n > 0) {
295: PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id));
296: PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
297: } else {
298: /* Create null filespace to match null memspace. */
299: PetscStackCallHDF5Return(filespace,H5Screate,(H5S_NULL));
300: }
302: ISGetIndices(is, &ind);
303: PetscStackCallHDF5(H5Dwrite,(dset_id, inttype, memspace, filespace, hdf5->dxpl_id, ind));
304: PetscStackCallHDF5(H5Fflush,(file_id, H5F_SCOPE_GLOBAL));
305: ISRestoreIndices(is, &ind);
307: /* Close/release resources */
308: PetscStackCallHDF5(H5Gclose,(group));
309: PetscStackCallHDF5(H5Sclose,(filespace));
310: PetscStackCallHDF5(H5Sclose,(memspace));
311: PetscStackCallHDF5(H5Dclose,(dset_id));
312: PetscInfo1(is, "Wrote IS object with name %s\n", isname);
313: return(0);
314: }
315: #endif
317: static PetscErrorCode ISView_General_Binary(IS is,PetscViewer viewer)
318: {
320: PetscBool skipHeader,useMPIIO;
321: IS_General *isa = (IS_General*) is->data;
322: PetscMPIInt rank,size,mesgsize,tag = ((PetscObject)viewer)->tag, mesglen;
323: PetscInt n,N,len,j,tr[2];
324: int fdes;
325: MPI_Status status;
326: PetscInt message_count,flowcontrolcount,*values;
329: /* ISGetLayout(is,&map); */
330: PetscLayoutGetLocalSize(is->map, &n);
331: PetscLayoutGetSize(is->map, &N);
333: tr[0] = IS_FILE_CLASSID;
334: tr[1] = N;
336: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
337: if (!skipHeader) {
338: PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_FALSE);
339: }
341: PetscViewerBinaryGetUseMPIIO(viewer,&useMPIIO);
342: #if defined(PETSC_HAVE_MPIIO)
343: if (useMPIIO) {
344: MPI_File mfdes;
345: MPI_Offset off;
346: PetscMPIInt lsize;
347: PetscInt rstart;
348: const PetscInt *iarray;
350: PetscMPIIntCast(n,&lsize);
351: PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);
352: PetscViewerBinaryGetMPIIOOffset(viewer,&off);
353: PetscLayoutGetRange(is->map,&rstart,NULL);
354: off += rstart*(MPI_Offset)sizeof(PetscInt); /* off is MPI_Offset, not PetscMPIInt */
355: MPI_File_set_view(mfdes,off,MPIU_INT,MPIU_INT,(char*)"native",MPI_INFO_NULL);
356: ISGetIndices(is,&iarray);
357: MPIU_File_write_all(mfdes,(void*)iarray,lsize,MPIU_INT,MPI_STATUS_IGNORE);
358: ISRestoreIndices(is,&iarray);
359: PetscViewerBinaryAddMPIIOOffset(viewer,N*(MPI_Offset)sizeof(PetscInt));
360: return(0);
361: }
362: #endif
364: PetscViewerBinaryGetDescriptor(viewer,&fdes);
365: MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
366: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
368: /* determine maximum message to arrive */
369: MPI_Reduce(&n,&len,1,MPIU_INT,MPI_MAX,0,PetscObjectComm((PetscObject)is));
371: PetscViewerFlowControlStart(viewer,&message_count,&flowcontrolcount);
372: if (!rank) {
373: PetscBinaryWrite(fdes,isa->idx,n,PETSC_INT,PETSC_FALSE);
375: PetscMalloc1(len,&values);
376: PetscMPIIntCast(len,&mesgsize);
377: /* receive and save messages */
378: for (j=1; j<size; j++) {
379: PetscViewerFlowControlStepMaster(viewer,j,&message_count,flowcontrolcount);
380: MPI_Recv(values,mesgsize,MPIU_INT,j,tag,PetscObjectComm((PetscObject)is),&status);
381: MPI_Get_count(&status,MPIU_INT,&mesglen);
382: PetscBinaryWrite(fdes,values,(PetscInt)mesglen,PETSC_INT,PETSC_TRUE);
383: }
384: PetscViewerFlowControlEndMaster(viewer,&message_count);
385: PetscFree(values);
386: } else {
387: PetscViewerFlowControlStepWorker(viewer,rank,&message_count);
388: PetscMPIIntCast(n,&mesgsize);
389: MPI_Send(isa->idx,mesgsize,MPIU_INT,0,tag,PetscObjectComm((PetscObject)is));
390: PetscViewerFlowControlEndWorker(viewer,&message_count);
391: }
392: return(0);
393: }
395: static PetscErrorCode ISView_General(IS is,PetscViewer viewer)
396: {
397: IS_General *sub = (IS_General*)is->data;
399: PetscInt i,n,*idx = sub->idx;
400: PetscBool iascii,isbinary,ishdf5;
403: PetscLayoutGetLocalSize(is->map, &n);
404: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
405: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
406: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
407: if (iascii) {
408: MPI_Comm comm;
409: PetscMPIInt rank,size;
410: PetscViewerFormat fmt;
411: PetscBool isperm;
413: PetscObjectGetComm((PetscObject)viewer,&comm);
414: MPI_Comm_rank(comm,&rank);
415: MPI_Comm_size(comm,&size);
417: PetscViewerGetFormat(viewer,&fmt);
418: ISGetInfo(is,IS_PERMUTATION,IS_LOCAL,PETSC_FALSE,&isperm);
419: if (isperm && fmt != PETSC_VIEWER_ASCII_MATLAB) {PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");}
420: PetscViewerASCIIPushSynchronized(viewer);
421: if (size > 1) {
422: if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
423: const char* name;
425: PetscObjectGetName((PetscObject)is,&name);
426: PetscViewerASCIISynchronizedPrintf(viewer,"%s_%d = [...\n",name,rank);
427: for (i=0; i<n; i++) {
428: PetscViewerASCIISynchronizedPrintf(viewer,"%D\n",idx[i]+1);
429: }
430: PetscViewerASCIISynchronizedPrintf(viewer,"];\n");
431: } else {
432: PetscInt st = 0;
434: if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
435: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in set %D\n",rank,n);
436: for (i=0; i<n; i++) {
437: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i + st,idx[i]);
438: }
439: }
440: } else {
441: if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
442: const char* name;
444: PetscObjectGetName((PetscObject)is,&name);
445: PetscViewerASCIISynchronizedPrintf(viewer,"%s = [...\n",name);
446: for (i=0; i<n; i++) {
447: PetscViewerASCIISynchronizedPrintf(viewer,"%D\n",idx[i]+1);
448: }
449: PetscViewerASCIISynchronizedPrintf(viewer,"];\n");
450: } else {
451: PetscInt st = 0;
453: if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
454: PetscViewerASCIISynchronizedPrintf(viewer,"Number of indices in set %D\n",n);
455: for (i=0; i<n; i++) {
456: PetscViewerASCIISynchronizedPrintf(viewer,"%D %D\n",i + st,idx[i]);
457: }
458: }
459: }
460: PetscViewerFlush(viewer);
461: PetscViewerASCIIPopSynchronized(viewer);
462: } else if (isbinary) {
463: ISView_General_Binary(is,viewer);
464: } else if (ishdf5) {
465: #if defined(PETSC_HAVE_HDF5)
466: ISView_General_HDF5(is,viewer);
467: #endif
468: }
469: return(0);
470: }
472: static PetscErrorCode ISSort_General(IS is)
473: {
474: IS_General *sub = (IS_General*)is->data;
475: PetscInt n;
479: PetscLayoutGetLocalSize(is->map, &n);
480: PetscSortInt(n,sub->idx);
481: return(0);
482: }
484: static PetscErrorCode ISSortRemoveDups_General(IS is)
485: {
486: IS_General *sub = (IS_General*)is->data;
487: PetscLayout map;
488: PetscInt n;
489: PetscBool sorted;
493: PetscLayoutGetLocalSize(is->map, &n);
494: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sorted);
495: if (sorted) {
496: PetscSortedRemoveDupsInt(&n,sub->idx);
497: } else {
498: PetscSortRemoveDupsInt(&n,sub->idx);
499: }
500: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map);
501: PetscLayoutDestroy(&is->map);
502: is->map = map;
503: return(0);
504: }
506: static PetscErrorCode ISSorted_General(IS is,PetscBool *flg)
507: {
511: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,flg);
512: return(0);
513: }
515: PetscErrorCode ISToGeneral_General(IS is)
516: {
518: return(0);
519: }
521: static struct _ISOps myops = { ISGetIndices_General,
522: ISRestoreIndices_General,
523: ISInvertPermutation_General,
524: ISSort_General,
525: ISSortRemoveDups_General,
526: ISSorted_General,
527: ISDuplicate_General,
528: ISDestroy_General,
529: ISView_General,
530: ISLoad_Default,
531: ISCopy_General,
532: ISToGeneral_General,
533: ISOnComm_General,
534: ISSetBlockSize_General,
535: ISContiguousLocal_General,
536: ISLocate_General,
537: /* no specializations of {sorted,unique,perm,interval}{local,global}
538: * because the default checks in ISGetInfo_XXX in index.c are exactly
539: * what we would do for ISGeneral */
540: NULL,
541: NULL,
542: NULL,
543: NULL,
544: NULL,
545: NULL,
546: NULL,
547: NULL};
549: PETSC_INTERN PetscErrorCode ISSetUp_General(IS);
551: PetscErrorCode ISSetUp_General(IS is)
552: {
554: IS_General *sub = (IS_General*)is->data;
555: const PetscInt *idx = sub->idx;
556: PetscInt n,i,min,max;
559: PetscLayoutGetLocalSize(is->map, &n);
561: if (n) {
562: min = max = idx[0];
563: for (i=1; i<n; i++) {
564: if (idx[i] < min) min = idx[i];
565: if (idx[i] > max) max = idx[i];
566: }
567: is->min = min;
568: is->max = max;
569: } else {
570: is->min = PETSC_MAX_INT;
571: is->max = PETSC_MIN_INT;
572: }
573: return(0);
574: }
576: /*@
577: ISCreateGeneral - Creates a data structure for an index set
578: containing a list of integers.
580: Collective
582: Input Parameters:
583: + comm - the MPI communicator
584: . n - the length of the index set
585: . idx - the list of integers
586: - mode - PETSC_COPY_VALUES, PETSC_OWN_POINTER, or PETSC_USE_POINTER; see PetscCopyMode for meaning of this flag.
588: Output Parameter:
589: . is - the new index set
591: Notes:
592: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
593: conceptually the same as MPI_Group operations. The IS are then
594: distributed sets of indices and thus certain operations on them are
595: collective.
598: Level: beginner
601: .seealso: ISCreateStride(), ISCreateBlock(), ISAllGather(), PETSC_COPY_VALUES, PETSC_OWN_POINTER, PETSC_USE_POINTER, PetscCopyMode
602: @*/
603: PetscErrorCode ISCreateGeneral(MPI_Comm comm,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
604: {
608: ISCreate(comm,is);
609: ISSetType(*is,ISGENERAL);
610: ISGeneralSetIndices(*is,n,idx,mode);
611: return(0);
612: }
614: /*@
615: ISGeneralSetIndices - Sets the indices for an ISGENERAL index set
617: Collective on IS
619: Input Parameters:
620: + is - the index set
621: . n - the length of the index set
622: . idx - the list of integers
623: - mode - see PetscCopyMode for meaning of this flag.
625: Level: beginner
628: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
629: @*/
630: PetscErrorCode ISGeneralSetIndices(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
631: {
635: ISClearInfoCache(is,PETSC_FALSE);
636: PetscUseMethod(is,"ISGeneralSetIndices_C",(IS,PetscInt,const PetscInt[],PetscCopyMode),(is,n,idx,mode));
637: return(0);
638: }
640: PetscErrorCode ISGeneralSetIndices_General(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
641: {
642: PetscLayout map;
644: IS_General *sub = (IS_General*)is->data;
647: if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
650: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is),n,PETSC_DECIDE,is->map->bs,&map);
651: PetscLayoutDestroy(&is->map);
652: is->map = map;
654: if (sub->allocated) {PetscFree(sub->idx);}
655: if (mode == PETSC_COPY_VALUES) {
656: PetscMalloc1(n,&sub->idx);
657: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
658: PetscArraycpy(sub->idx,idx,n);
659: sub->allocated = PETSC_TRUE;
660: } else if (mode == PETSC_OWN_POINTER) {
661: sub->idx = (PetscInt*)idx;
662: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
663: sub->allocated = PETSC_TRUE;
664: } else {
665: sub->idx = (PetscInt*)idx;
666: sub->allocated = PETSC_FALSE;
667: }
669: ISSetUp_General(is);
670: ISViewFromOptions(is,NULL,"-is_view");
671: return(0);
672: }
674: static PetscErrorCode ISGeneralFilter_General(IS is, PetscInt start, PetscInt end)
675: {
676: IS_General *sub = (IS_General*)is->data;
677: PetscInt *idx = sub->idx,*idxnew;
678: PetscInt i,n = is->map->n,nnew = 0,o;
682: for (i=0; i<n; ++i)
683: if (idx[i] >= start && idx[i] < end)
684: nnew++;
685: PetscMalloc1(nnew, &idxnew);
686: for (o=0, i=0; i<n; i++) {
687: if (idx[i] >= start && idx[i] < end)
688: idxnew[o++] = idx[i];
689: }
690: ISGeneralSetIndices_General(is,nnew,idxnew,PETSC_OWN_POINTER);
691: return(0);
692: }
694: /*@
695: ISGeneralFilter - Remove all points outside of [start, end)
697: Collective on IS
699: Input Parameters:
700: + is - the index set
701: . start - the lowest index kept
702: - end - one more than the highest index kept
704: Level: beginner
706: .seealso: ISCreateGeneral(), ISGeneralSetIndices()
707: @*/
708: PetscErrorCode ISGeneralFilter(IS is, PetscInt start, PetscInt end)
709: {
714: ISClearInfoCache(is,PETSC_FALSE);
715: PetscUseMethod(is,"ISGeneralFilter_C",(IS,PetscInt,PetscInt),(is,start,end));
716: return(0);
717: }
719: PETSC_EXTERN PetscErrorCode ISCreate_General(IS is)
720: {
722: IS_General *sub;
725: PetscNewLog(is,&sub);
726: is->data = (void *) sub;
727: PetscMemcpy(is->ops,&myops,sizeof(myops));
728: PetscObjectComposeFunction((PetscObject)is,"ISGeneralSetIndices_C",ISGeneralSetIndices_General);
729: PetscObjectComposeFunction((PetscObject)is,"ISGeneralFilter_C",ISGeneralFilter_General);
730: return(0);
731: }