Actual source code: general.c
petsc-main 2021-04-20
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: /* F90Array1dCreate() inside ISRestoreArrayF90() does not keep array when zero length array */
137: if (in->map->n > 0 && *idx != sub->idx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
138: return(0);
139: }
141: static PetscErrorCode ISInvertPermutation_General(IS is,PetscInt nlocal,IS *isout)
142: {
143: IS_General *sub = (IS_General*)is->data;
144: PetscInt i,*ii,n,nstart;
145: const PetscInt *idx = sub->idx;
146: PetscMPIInt size;
147: IS istmp,nistmp;
151: PetscLayoutGetLocalSize(is->map, &n);
152: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
153: if (size == 1) {
154: PetscMalloc1(n,&ii);
155: for (i=0; i<n; i++) ii[idx[i]] = i;
156: ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_OWN_POINTER,isout);
157: ISSetPermutation(*isout);
158: } else {
159: /* crude, nonscalable get entire IS on each processor */
160: ISAllGather(is,&istmp);
161: ISSetPermutation(istmp);
162: ISInvertPermutation(istmp,PETSC_DECIDE,&nistmp);
163: ISDestroy(&istmp);
164: /* get the part we need */
165: if (nlocal == PETSC_DECIDE) nlocal = n;
166: MPI_Scan(&nlocal,&nstart,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)is));
167: if (PetscDefined(USE_DEBUG)) {
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: nstart -= nlocal;
177: ISGetIndices(nistmp,&idx);
178: ISCreateGeneral(PetscObjectComm((PetscObject)is),nlocal,idx+nstart,PETSC_COPY_VALUES,isout);
179: ISRestoreIndices(nistmp,&idx);
180: ISDestroy(&nistmp);
181: }
182: return(0);
183: }
185: #if defined(PETSC_HAVE_HDF5)
186: static PetscErrorCode ISView_General_HDF5(IS is, PetscViewer viewer)
187: {
188: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
189: hid_t filespace; /* file dataspace identifier */
190: hid_t chunkspace; /* chunk dataset property identifier */
191: hid_t dset_id; /* dataset identifier */
192: hid_t memspace; /* memory dataspace identifier */
193: hid_t inttype; /* int type (H5T_NATIVE_INT or H5T_NATIVE_LLONG) */
194: hid_t file_id, group;
195: hsize_t dim, maxDims[3], dims[3], chunkDims[3], count[3],offset[3];
196: PetscInt bs, N, n, timestep, low;
197: const PetscInt *ind;
198: const char *isname;
199: PetscErrorCode ierr;
202: ISGetBlockSize(is,&bs);
203: bs = PetscMax(bs, 1); /* If N = 0, bs = 0 as well */
204: PetscViewerHDF5OpenGroup(viewer, &file_id, &group);
205: PetscViewerHDF5GetTimestep(viewer, ×tep);
207: /* Create the dataspace for the dataset.
208: *
209: * dims - holds the current dimensions of the dataset
210: *
211: * maxDims - holds the maximum dimensions of the dataset (unlimited
212: * for the number of time steps with the current dimensions for the
213: * other dimensions; so only additional time steps can be added).
214: *
215: * chunkDims - holds the size of a single time step (required to
216: * permit extending dataset).
217: */
218: dim = 0;
219: if (timestep >= 0) {
220: dims[dim] = timestep+1;
221: maxDims[dim] = H5S_UNLIMITED;
222: chunkDims[dim] = 1;
223: ++dim;
224: }
225: ISGetSize(is, &N);
226: ISGetLocalSize(is, &n);
227: PetscHDF5IntCast(N/bs,dims + dim);
229: maxDims[dim] = dims[dim];
230: chunkDims[dim] = PetscMax(1,dims[dim]);
231: ++dim;
232: if (bs >= 1) {
233: dims[dim] = bs;
234: maxDims[dim] = dims[dim];
235: chunkDims[dim] = dims[dim];
236: ++dim;
237: }
238: PetscStackCallHDF5Return(filespace,H5Screate_simple,(dim, dims, maxDims));
240: #if defined(PETSC_USE_64BIT_INDICES)
241: inttype = H5T_NATIVE_LLONG;
242: #else
243: inttype = H5T_NATIVE_INT;
244: #endif
246: /* Create the dataset with default properties and close filespace */
247: PetscObjectGetName((PetscObject) is, &isname);
248: if (!H5Lexists(group, isname, H5P_DEFAULT)) {
249: /* Create chunk */
250: PetscStackCallHDF5Return(chunkspace,H5Pcreate,(H5P_DATASET_CREATE));
251: PetscStackCallHDF5(H5Pset_chunk,(chunkspace, dim, chunkDims));
253: PetscStackCallHDF5Return(dset_id,H5Dcreate2,(group, isname, inttype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
254: PetscStackCallHDF5(H5Pclose,(chunkspace));
255: } else {
256: PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, isname, H5P_DEFAULT));
257: PetscStackCallHDF5(H5Dset_extent,(dset_id, dims));
258: }
259: PetscStackCallHDF5(H5Sclose,(filespace));
261: /* Each process defines a dataset and writes it to the hyperslab in the file */
262: dim = 0;
263: if (timestep >= 0) {
264: count[dim] = 1;
265: ++dim;
266: }
267: PetscHDF5IntCast(n/bs,count + dim);
268: ++dim;
269: if (bs >= 1) {
270: count[dim] = bs;
271: ++dim;
272: }
273: if (n > 0 || H5_VERSION_GE(1,10,0)) {
274: PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL));
275: } else {
276: /* Can't create dataspace with zero for any dimension, so create null dataspace. */
277: PetscStackCallHDF5Return(memspace,H5Screate,(H5S_NULL));
278: }
280: /* Select hyperslab in the file */
281: PetscLayoutGetRange(is->map, &low, NULL);
282: dim = 0;
283: if (timestep >= 0) {
284: offset[dim] = timestep;
285: ++dim;
286: }
287: PetscHDF5IntCast(low/bs,offset + dim);
288: ++dim;
289: if (bs >= 1) {
290: offset[dim] = 0;
291: ++dim;
292: }
293: if (n > 0 || H5_VERSION_GE(1,10,0)) {
294: PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id));
295: PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
296: } else {
297: /* Create null filespace to match null memspace. */
298: PetscStackCallHDF5Return(filespace,H5Screate,(H5S_NULL));
299: }
301: ISGetIndices(is, &ind);
302: PetscStackCallHDF5(H5Dwrite,(dset_id, inttype, memspace, filespace, hdf5->dxpl_id, ind));
303: PetscStackCallHDF5(H5Fflush,(file_id, H5F_SCOPE_GLOBAL));
304: ISRestoreIndices(is, &ind);
306: /* Close/release resources */
307: PetscStackCallHDF5(H5Gclose,(group));
308: PetscStackCallHDF5(H5Sclose,(filespace));
309: PetscStackCallHDF5(H5Sclose,(memspace));
310: PetscStackCallHDF5(H5Dclose,(dset_id));
311: PetscInfo1(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;
320: PetscInt i,n,*idx = sub->idx;
321: PetscBool iascii,isbinary,ishdf5;
324: PetscLayoutGetLocalSize(is->map, &n);
325: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
326: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
327: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
328: if (iascii) {
329: MPI_Comm comm;
330: PetscMPIInt rank,size;
331: PetscViewerFormat fmt;
332: PetscBool isperm;
334: PetscObjectGetComm((PetscObject)viewer,&comm);
335: MPI_Comm_rank(comm,&rank);
336: MPI_Comm_size(comm,&size);
338: PetscViewerGetFormat(viewer,&fmt);
339: ISGetInfo(is,IS_PERMUTATION,IS_LOCAL,PETSC_FALSE,&isperm);
340: if (isperm && fmt != PETSC_VIEWER_ASCII_MATLAB) {PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");}
341: PetscViewerASCIIPushSynchronized(viewer);
342: if (size > 1) {
343: if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
344: const char* name;
346: PetscObjectGetName((PetscObject)is,&name);
347: PetscViewerASCIISynchronizedPrintf(viewer,"%s_%d = [...\n",name,rank);
348: for (i=0; i<n; i++) {
349: PetscViewerASCIISynchronizedPrintf(viewer,"%D\n",idx[i]+1);
350: }
351: PetscViewerASCIISynchronizedPrintf(viewer,"];\n");
352: } else {
353: PetscInt st = 0;
355: if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
356: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in set %D\n",rank,n);
357: for (i=0; i<n; i++) {
358: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i + st,idx[i]);
359: }
360: }
361: } else {
362: if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
363: const char* name;
365: PetscObjectGetName((PetscObject)is,&name);
366: PetscViewerASCIISynchronizedPrintf(viewer,"%s = [...\n",name);
367: for (i=0; i<n; i++) {
368: PetscViewerASCIISynchronizedPrintf(viewer,"%D\n",idx[i]+1);
369: }
370: PetscViewerASCIISynchronizedPrintf(viewer,"];\n");
371: } else {
372: PetscInt st = 0;
374: if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
375: PetscViewerASCIISynchronizedPrintf(viewer,"Number of indices in set %D\n",n);
376: for (i=0; i<n; i++) {
377: PetscViewerASCIISynchronizedPrintf(viewer,"%D %D\n",i + st,idx[i]);
378: }
379: }
380: }
381: PetscViewerFlush(viewer);
382: PetscViewerASCIIPopSynchronized(viewer);
383: } else if (isbinary) {
384: ISView_Binary(is,viewer);
385: } else if (ishdf5) {
386: #if defined(PETSC_HAVE_HDF5)
387: ISView_General_HDF5(is,viewer);
388: #endif
389: }
390: return(0);
391: }
393: static PetscErrorCode ISSort_General(IS is)
394: {
395: IS_General *sub = (IS_General*)is->data;
396: PetscInt n;
400: PetscLayoutGetLocalSize(is->map, &n);
401: PetscIntSortSemiOrdered(n,sub->idx);
402: return(0);
403: }
405: static PetscErrorCode ISSortRemoveDups_General(IS is)
406: {
407: IS_General *sub = (IS_General*)is->data;
408: PetscLayout map;
409: PetscInt n;
410: PetscBool sorted;
414: PetscLayoutGetLocalSize(is->map, &n);
415: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sorted);
416: if (sorted) {
417: PetscSortedRemoveDupsInt(&n,sub->idx);
418: } else {
419: PetscSortRemoveDupsInt(&n,sub->idx);
420: }
421: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map);
422: PetscLayoutDestroy(&is->map);
423: is->map = map;
424: return(0);
425: }
427: static PetscErrorCode ISSorted_General(IS is,PetscBool *flg)
428: {
432: ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,flg);
433: return(0);
434: }
436: PetscErrorCode ISToGeneral_General(IS is)
437: {
439: return(0);
440: }
442: static struct _ISOps myops = { ISGetIndices_General,
443: ISRestoreIndices_General,
444: ISInvertPermutation_General,
445: ISSort_General,
446: ISSortRemoveDups_General,
447: ISSorted_General,
448: ISDuplicate_General,
449: ISDestroy_General,
450: ISView_General,
451: ISLoad_Default,
452: ISCopy_General,
453: ISToGeneral_General,
454: ISOnComm_General,
455: ISSetBlockSize_General,
456: ISContiguousLocal_General,
457: ISLocate_General,
458: /* no specializations of {sorted,unique,perm,interval}{local,global}
459: * because the default checks in ISGetInfo_XXX in index.c are exactly
460: * what we would do for ISGeneral */
461: NULL,
462: NULL,
463: NULL,
464: NULL,
465: NULL,
466: NULL,
467: NULL,
468: NULL};
470: PETSC_INTERN PetscErrorCode ISSetUp_General(IS);
472: PetscErrorCode ISSetUp_General(IS is)
473: {
475: IS_General *sub = (IS_General*)is->data;
476: const PetscInt *idx = sub->idx;
477: PetscInt n,i,min,max;
480: PetscLayoutGetLocalSize(is->map, &n);
482: if (n) {
483: min = max = idx[0];
484: for (i=1; i<n; i++) {
485: if (idx[i] < min) min = idx[i];
486: if (idx[i] > max) max = idx[i];
487: }
488: is->min = min;
489: is->max = max;
490: } else {
491: is->min = PETSC_MAX_INT;
492: is->max = PETSC_MIN_INT;
493: }
494: return(0);
495: }
497: /*@
498: ISCreateGeneral - Creates a data structure for an index set
499: containing a list of integers.
501: Collective
503: Input Parameters:
504: + comm - the MPI communicator
505: . n - the length of the index set
506: . idx - the list of integers
507: - mode - PETSC_COPY_VALUES, PETSC_OWN_POINTER, or PETSC_USE_POINTER; see PetscCopyMode for meaning of this flag.
509: Output Parameter:
510: . is - the new index set
512: Notes:
513: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
514: conceptually the same as MPI_Group operations. The IS are then
515: distributed sets of indices and thus certain operations on them are
516: collective.
518: Level: beginner
520: .seealso: ISCreateStride(), ISCreateBlock(), ISAllGather(), PETSC_COPY_VALUES, PETSC_OWN_POINTER, PETSC_USE_POINTER, PetscCopyMode
521: @*/
522: PetscErrorCode ISCreateGeneral(MPI_Comm comm,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
523: {
527: ISCreate(comm,is);
528: ISSetType(*is,ISGENERAL);
529: ISGeneralSetIndices(*is,n,idx,mode);
530: return(0);
531: }
533: /*@
534: ISGeneralSetIndices - Sets the indices for an ISGENERAL index set
536: Collective on IS
538: Input Parameters:
539: + is - the index set
540: . n - the length of the index set
541: . idx - the list of integers
542: - mode - see PetscCopyMode for meaning of this flag.
544: Level: beginner
546: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather(), ISBlockSetIndices(), ISGENERAL, PetscCopyMode
547: @*/
548: PetscErrorCode ISGeneralSetIndices(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
549: {
553: ISClearInfoCache(is,PETSC_FALSE);
554: PetscUseMethod(is,"ISGeneralSetIndices_C",(IS,PetscInt,const PetscInt[],PetscCopyMode),(is,n,idx,mode));
555: return(0);
556: }
558: PetscErrorCode ISGeneralSetIndices_General(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
559: {
560: PetscLayout map;
562: IS_General *sub = (IS_General*)is->data;
565: if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
568: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is),n,PETSC_DECIDE,is->map->bs,&map);
569: PetscLayoutDestroy(&is->map);
570: is->map = map;
572: if (sub->allocated) {PetscFree(sub->idx);}
573: if (mode == PETSC_COPY_VALUES) {
574: PetscMalloc1(n,&sub->idx);
575: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
576: PetscArraycpy(sub->idx,idx,n);
577: sub->allocated = PETSC_TRUE;
578: } else if (mode == PETSC_OWN_POINTER) {
579: sub->idx = (PetscInt*)idx;
580: PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
581: sub->allocated = PETSC_TRUE;
582: } else {
583: sub->idx = (PetscInt*)idx;
584: sub->allocated = PETSC_FALSE;
585: }
587: ISSetUp_General(is);
588: ISViewFromOptions(is,NULL,"-is_view");
589: return(0);
590: }
592: static PetscErrorCode ISGeneralFilter_General(IS is, PetscInt start, PetscInt end)
593: {
594: IS_General *sub = (IS_General*)is->data;
595: PetscInt *idx = sub->idx,*idxnew;
596: PetscInt i,n = is->map->n,nnew = 0,o;
600: for (i=0; i<n; ++i)
601: if (idx[i] >= start && idx[i] < end)
602: nnew++;
603: PetscMalloc1(nnew, &idxnew);
604: for (o=0, i=0; i<n; i++) {
605: if (idx[i] >= start && idx[i] < end)
606: idxnew[o++] = idx[i];
607: }
608: ISGeneralSetIndices_General(is,nnew,idxnew,PETSC_OWN_POINTER);
609: return(0);
610: }
612: /*@
613: ISGeneralFilter - Remove all indices outside of [start, end)
615: Collective on IS
617: Input Parameters:
618: + is - the index set
619: . start - the lowest index kept
620: - end - one more than the highest index kept
622: Level: beginner
624: .seealso: ISCreateGeneral(), ISGeneralSetIndices()
625: @*/
626: PetscErrorCode ISGeneralFilter(IS is, PetscInt start, PetscInt end)
627: {
632: ISClearInfoCache(is,PETSC_FALSE);
633: PetscUseMethod(is,"ISGeneralFilter_C",(IS,PetscInt,PetscInt),(is,start,end));
634: return(0);
635: }
637: PETSC_EXTERN PetscErrorCode ISCreate_General(IS is)
638: {
640: IS_General *sub;
643: PetscNewLog(is,&sub);
644: is->data = (void *) sub;
645: PetscMemcpy(is->ops,&myops,sizeof(myops));
646: PetscObjectComposeFunction((PetscObject)is,"ISGeneralSetIndices_C",ISGeneralSetIndices_General);
647: PetscObjectComposeFunction((PetscObject)is,"ISGeneralFilter_C",ISGeneralFilter_General);
648: return(0);
649: }