Actual source code: general.c
petsc-master 2020-08-06
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 (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: PetscInt bs, N, n, timestep, low;
196: const PetscInt *ind;
197: const char *isname;
198: PetscErrorCode ierr;
201: ISGetBlockSize(is,&bs);
202: bs = PetscMax(bs, 1); /* If N = 0, bs = 0 as well */
203: PetscViewerHDF5OpenGroup(viewer, &file_id, &group);
204: PetscViewerHDF5GetTimestep(viewer, ×tep);
206: /* Create the dataspace for the dataset.
207: *
208: * dims - holds the current dimensions of the dataset
209: *
210: * maxDims - holds the maximum dimensions of the dataset (unlimited
211: * for the number of time steps with the current dimensions for the
212: * other dimensions; so only additional time steps can be added).
213: *
214: * chunkDims - holds the size of a single time step (required to
215: * permit extending dataset).
216: */
217: dim = 0;
218: if (timestep >= 0) {
219: dims[dim] = timestep+1;
220: maxDims[dim] = H5S_UNLIMITED;
221: chunkDims[dim] = 1;
222: ++dim;
223: }
224: ISGetSize(is, &N);
225: ISGetLocalSize(is, &n);
226: PetscHDF5IntCast(N/bs,dims + dim);
228: maxDims[dim] = dims[dim];
229: chunkDims[dim] = PetscMax(1,dims[dim]);
230: ++dim;
231: if (bs >= 1) {
232: dims[dim] = bs;
233: maxDims[dim] = dims[dim];
234: chunkDims[dim] = dims[dim];
235: ++dim;
236: }
237: PetscStackCallHDF5Return(filespace,H5Screate_simple,(dim, dims, maxDims));
239: #if defined(PETSC_USE_64BIT_INDICES)
240: inttype = H5T_NATIVE_LLONG;
241: #else
242: inttype = H5T_NATIVE_INT;
243: #endif
245: /* Create the dataset with default properties and close filespace */
246: PetscObjectGetName((PetscObject) is, &isname);
247: if (!H5Lexists(group, isname, H5P_DEFAULT)) {
248: /* Create chunk */
249: PetscStackCallHDF5Return(chunkspace,H5Pcreate,(H5P_DATASET_CREATE));
250: PetscStackCallHDF5(H5Pset_chunk,(chunkspace, dim, chunkDims));
252: PetscStackCallHDF5Return(dset_id,H5Dcreate2,(group, isname, inttype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
253: PetscStackCallHDF5(H5Pclose,(chunkspace));
254: } else {
255: PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, isname, H5P_DEFAULT));
256: PetscStackCallHDF5(H5Dset_extent,(dset_id, dims));
257: }
258: PetscStackCallHDF5(H5Sclose,(filespace));
260: /* Each process defines a dataset and writes it to the hyperslab in the file */
261: dim = 0;
262: if (timestep >= 0) {
263: count[dim] = 1;
264: ++dim;
265: }
266: PetscHDF5IntCast(n/bs,count + dim);
267: ++dim;
268: if (bs >= 1) {
269: count[dim] = bs;
270: ++dim;
271: }
272: if (n > 0 || H5_VERSION_GE(1,10,0)) {
273: PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL));
274: } else {
275: /* Can't create dataspace with zero for any dimension, so create null dataspace. */
276: PetscStackCallHDF5Return(memspace,H5Screate,(H5S_NULL));
277: }
279: /* Select hyperslab in the file */
280: PetscLayoutGetRange(is->map, &low, NULL);
281: dim = 0;
282: if (timestep >= 0) {
283: offset[dim] = timestep;
284: ++dim;
285: }
286: PetscHDF5IntCast(low/bs,offset + dim);
287: ++dim;
288: if (bs >= 1) {
289: offset[dim] = 0;
290: ++dim;
291: }
292: if (n > 0 || H5_VERSION_GE(1,10,0)) {
293: PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id));
294: PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
295: } else {
296: /* Create null filespace to match null memspace. */
297: PetscStackCallHDF5Return(filespace,H5Screate,(H5S_NULL));
298: }
300: ISGetIndices(is, &ind);
301: PetscStackCallHDF5(H5Dwrite,(dset_id, inttype, memspace, filespace, hdf5->dxpl_id, ind));
302: PetscStackCallHDF5(H5Fflush,(file_id, H5F_SCOPE_GLOBAL));
303: ISRestoreIndices(is, &ind);
305: /* Close/release resources */
306: PetscStackCallHDF5(H5Gclose,(group));
307: PetscStackCallHDF5(H5Sclose,(filespace));
308: PetscStackCallHDF5(H5Sclose,(memspace));
309: PetscStackCallHDF5(H5Dclose,(dset_id));
310: PetscInfo1(is, "Wrote IS object with name %s\n", isname);
311: return(0);
312: }
313: #endif
315: static PetscErrorCode ISView_General(IS is,PetscViewer viewer)
316: {
317: IS_General *sub = (IS_General*)is->data;
319: PetscInt i,n,*idx = sub->idx;
320: PetscBool iascii,isbinary,ishdf5;
323: PetscLayoutGetLocalSize(is->map, &n);
324: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
325: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
326: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
327: if (iascii) {
328: MPI_Comm comm;
329: PetscMPIInt rank,size;
330: PetscViewerFormat fmt;
331: PetscBool isperm;
333: PetscObjectGetComm((PetscObject)viewer,&comm);
334: MPI_Comm_rank(comm,&rank);
335: MPI_Comm_size(comm,&size);
337: PetscViewerGetFormat(viewer,&fmt);
338: ISGetInfo(is,IS_PERMUTATION,IS_LOCAL,PETSC_FALSE,&isperm);
339: if (isperm && fmt != PETSC_VIEWER_ASCII_MATLAB) {PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");}
340: PetscViewerASCIIPushSynchronized(viewer);
341: if (size > 1) {
342: if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
343: const char* name;
345: PetscObjectGetName((PetscObject)is,&name);
346: PetscViewerASCIISynchronizedPrintf(viewer,"%s_%d = [...\n",name,rank);
347: for (i=0; i<n; i++) {
348: PetscViewerASCIISynchronizedPrintf(viewer,"%D\n",idx[i]+1);
349: }
350: PetscViewerASCIISynchronizedPrintf(viewer,"];\n");
351: } else {
352: PetscInt st = 0;
354: if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
355: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in set %D\n",rank,n);
356: for (i=0; i<n; i++) {
357: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i + st,idx[i]);
358: }
359: }
360: } else {
361: if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
362: const char* name;
364: PetscObjectGetName((PetscObject)is,&name);
365: PetscViewerASCIISynchronizedPrintf(viewer,"%s = [...\n",name);
366: for (i=0; i<n; i++) {
367: PetscViewerASCIISynchronizedPrintf(viewer,"%D\n",idx[i]+1);
368: }
369: PetscViewerASCIISynchronizedPrintf(viewer,"];\n");
370: } else {
371: PetscInt st = 0;
373: if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
374: PetscViewerASCIISynchronizedPrintf(viewer,"Number of indices in set %D\n",n);
375: for (i=0; i<n; i++) {
376: PetscViewerASCIISynchronizedPrintf(viewer,"%D %D\n",i + st,idx[i]);
377: }
378: }
379: }
380: PetscViewerFlush(viewer);
381: PetscViewerASCIIPopSynchronized(viewer);
382: } else if (isbinary) {
383: ISView_Binary(is,viewer);
384: } else if (ishdf5) {
385: #if defined(PETSC_HAVE_HDF5)
386: ISView_General_HDF5(is,viewer);
387: #endif
388: }
389: return(0);
390: }
392: static PetscErrorCode ISSort_General(IS is)
393: {
394: IS_General *sub = (IS_General*)is->data;
395: PetscInt n;
399: PetscLayoutGetLocalSize(is->map, &n);
400: PetscSortInt(n,sub->idx);
401: return(0);
402: }
404: static PetscErrorCode ISSortRemoveDups_General(IS is)
405: {
406: IS_General *sub = (IS_General*)is->data;
407: PetscLayout map;
408: PetscInt n;
409: PetscBool sorted;
413: PetscLayoutGetLocalSize(is->map, &n);
414: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sorted);
415: if (sorted) {
416: PetscSortedRemoveDupsInt(&n,sub->idx);
417: } else {
418: PetscSortRemoveDupsInt(&n,sub->idx);
419: }
420: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map);
421: PetscLayoutDestroy(&is->map);
422: is->map = map;
423: return(0);
424: }
426: static PetscErrorCode ISSorted_General(IS is,PetscBool *flg)
427: {
431: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,flg);
432: return(0);
433: }
435: PetscErrorCode ISToGeneral_General(IS is)
436: {
438: return(0);
439: }
441: static struct _ISOps myops = { ISGetIndices_General,
442: ISRestoreIndices_General,
443: ISInvertPermutation_General,
444: ISSort_General,
445: ISSortRemoveDups_General,
446: ISSorted_General,
447: ISDuplicate_General,
448: ISDestroy_General,
449: ISView_General,
450: ISLoad_Default,
451: ISCopy_General,
452: ISToGeneral_General,
453: ISOnComm_General,
454: ISSetBlockSize_General,
455: ISContiguousLocal_General,
456: ISLocate_General,
457: /* no specializations of {sorted,unique,perm,interval}{local,global}
458: * because the default checks in ISGetInfo_XXX in index.c are exactly
459: * what we would do for ISGeneral */
460: NULL,
461: NULL,
462: NULL,
463: NULL,
464: NULL,
465: NULL,
466: NULL,
467: NULL};
469: PETSC_INTERN PetscErrorCode ISSetUp_General(IS);
471: PetscErrorCode ISSetUp_General(IS is)
472: {
474: IS_General *sub = (IS_General*)is->data;
475: const PetscInt *idx = sub->idx;
476: PetscInt n,i,min,max;
479: PetscLayoutGetLocalSize(is->map, &n);
481: if (n) {
482: min = max = idx[0];
483: for (i=1; i<n; i++) {
484: if (idx[i] < min) min = idx[i];
485: if (idx[i] > max) max = idx[i];
486: }
487: is->min = min;
488: is->max = max;
489: } else {
490: is->min = PETSC_MAX_INT;
491: is->max = PETSC_MIN_INT;
492: }
493: return(0);
494: }
496: /*@
497: ISCreateGeneral - Creates a data structure for an index set
498: containing a list of integers.
500: Collective
502: Input Parameters:
503: + comm - the MPI communicator
504: . n - the length of the index set
505: . idx - the list of integers
506: - mode - PETSC_COPY_VALUES, PETSC_OWN_POINTER, or PETSC_USE_POINTER; see PetscCopyMode for meaning of this flag.
508: Output Parameter:
509: . is - the new index set
511: Notes:
512: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
513: conceptually the same as MPI_Group operations. The IS are then
514: distributed sets of indices and thus certain operations on them are
515: collective.
518: Level: beginner
521: .seealso: ISCreateStride(), ISCreateBlock(), ISAllGather(), PETSC_COPY_VALUES, PETSC_OWN_POINTER, PETSC_USE_POINTER, PetscCopyMode
522: @*/
523: PetscErrorCode ISCreateGeneral(MPI_Comm comm,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
524: {
528: ISCreate(comm,is);
529: ISSetType(*is,ISGENERAL);
530: ISGeneralSetIndices(*is,n,idx,mode);
531: return(0);
532: }
534: /*@
535: ISGeneralSetIndices - Sets the indices for an ISGENERAL index set
537: Collective on IS
539: Input Parameters:
540: + is - the index set
541: . n - the length of the index set
542: . idx - the list of integers
543: - mode - see PetscCopyMode for meaning of this flag.
545: Level: beginner
548: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
549: @*/
550: PetscErrorCode ISGeneralSetIndices(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
551: {
555: ISClearInfoCache(is,PETSC_FALSE);
556: PetscUseMethod(is,"ISGeneralSetIndices_C",(IS,PetscInt,const PetscInt[],PetscCopyMode),(is,n,idx,mode));
557: return(0);
558: }
560: PetscErrorCode ISGeneralSetIndices_General(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
561: {
562: PetscLayout map;
564: IS_General *sub = (IS_General*)is->data;
567: if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
570: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is),n,PETSC_DECIDE,is->map->bs,&map);
571: PetscLayoutDestroy(&is->map);
572: is->map = map;
574: if (sub->allocated) {PetscFree(sub->idx);}
575: if (mode == PETSC_COPY_VALUES) {
576: PetscMalloc1(n,&sub->idx);
577: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
578: PetscArraycpy(sub->idx,idx,n);
579: sub->allocated = PETSC_TRUE;
580: } else if (mode == PETSC_OWN_POINTER) {
581: sub->idx = (PetscInt*)idx;
582: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
583: sub->allocated = PETSC_TRUE;
584: } else {
585: sub->idx = (PetscInt*)idx;
586: sub->allocated = PETSC_FALSE;
587: }
589: ISSetUp_General(is);
590: ISViewFromOptions(is,NULL,"-is_view");
591: return(0);
592: }
594: static PetscErrorCode ISGeneralFilter_General(IS is, PetscInt start, PetscInt end)
595: {
596: IS_General *sub = (IS_General*)is->data;
597: PetscInt *idx = sub->idx,*idxnew;
598: PetscInt i,n = is->map->n,nnew = 0,o;
602: for (i=0; i<n; ++i)
603: if (idx[i] >= start && idx[i] < end)
604: nnew++;
605: PetscMalloc1(nnew, &idxnew);
606: for (o=0, i=0; i<n; i++) {
607: if (idx[i] >= start && idx[i] < end)
608: idxnew[o++] = idx[i];
609: }
610: ISGeneralSetIndices_General(is,nnew,idxnew,PETSC_OWN_POINTER);
611: return(0);
612: }
614: /*@
615: ISGeneralFilter - Remove all points outside of [start, end)
617: Collective on IS
619: Input Parameters:
620: + is - the index set
621: . start - the lowest index kept
622: - end - one more than the highest index kept
624: Level: beginner
626: .seealso: ISCreateGeneral(), ISGeneralSetIndices()
627: @*/
628: PetscErrorCode ISGeneralFilter(IS is, PetscInt start, PetscInt end)
629: {
634: ISClearInfoCache(is,PETSC_FALSE);
635: PetscUseMethod(is,"ISGeneralFilter_C",(IS,PetscInt,PetscInt),(is,start,end));
636: return(0);
637: }
639: PETSC_EXTERN PetscErrorCode ISCreate_General(IS is)
640: {
642: IS_General *sub;
645: PetscNewLog(is,&sub);
646: is->data = (void *) sub;
647: PetscMemcpy(is->ops,&myops,sizeof(myops));
648: PetscObjectComposeFunction((PetscObject)is,"ISGeneralSetIndices_C",ISGeneralSetIndices_General);
649: PetscObjectComposeFunction((PetscObject)is,"ISGeneralFilter_C",ISGeneralFilter_General);
650: return(0);
651: }