Actual source code: general.c
petsc-master 2019-09-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 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, ×tep);
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: }