Actual source code: stride.c
petsc-3.14.6 2021-03-30
2: /*
3: Index sets of evenly space integers, defined by a
4: start, stride and length.
5: */
6: #include <petsc/private/isimpl.h>
7: #include <petscviewer.h>
9: typedef struct {
10: PetscInt first,step;
11: } IS_Stride;
13: static PetscErrorCode ISCopy_Stride(IS is,IS isy)
14: {
15: IS_Stride *is_stride = (IS_Stride*)is->data,*isy_stride = (IS_Stride*)isy->data;
19: PetscMemcpy(isy_stride,is_stride,sizeof(IS_Stride));
20: return(0);
21: }
23: PetscErrorCode ISDuplicate_Stride(IS is,IS *newIS)
24: {
26: IS_Stride *sub = (IS_Stride*)is->data;
29: ISCreateStride(PetscObjectComm((PetscObject)is),is->map->n,sub->first,sub->step,newIS);
30: return(0);
31: }
33: PetscErrorCode ISInvertPermutation_Stride(IS is,PetscInt nlocal,IS *perm)
34: {
35: PetscBool isident;
39: ISGetInfo(is,IS_IDENTITY,IS_GLOBAL,PETSC_TRUE,&isident);
40: if (isident) {
41: PetscInt rStart, rEnd;
43: PetscLayoutGetRange(is->map, &rStart, &rEnd);
44: ISCreateStride(PETSC_COMM_SELF,PetscMax(rEnd - rStart, 0),rStart,1,perm);
45: } else {
46: IS tmp;
47: const PetscInt *indices,n = is->map->n;
49: ISGetIndices(is,&indices);
50: ISCreateGeneral(PetscObjectComm((PetscObject)is),n,indices,PETSC_COPY_VALUES,&tmp);
51: ISSetPermutation(tmp);
52: ISRestoreIndices(is,&indices);
53: ISInvertPermutation(tmp,nlocal,perm);
54: ISDestroy(&tmp);
55: }
56: return(0);
57: }
59: /*@
60: ISStrideGetInfo - Returns the first index in a stride index set and the stride width.
62: Not Collective
64: Input Parameter:
65: . is - the index set
67: Output Parameters:
68: + first - the first index
69: - step - the stride width
71: Level: intermediate
73: Notes:
74: Returns info on stride index set. This is a pseudo-public function that
75: should not be needed by most users.
78: .seealso: ISCreateStride(), ISGetSize(), ISSTRIDE
79: @*/
80: PetscErrorCode ISStrideGetInfo(IS is,PetscInt *first,PetscInt *step)
81: {
82: IS_Stride *sub;
83: PetscBool flg;
90: PetscObjectTypeCompare((PetscObject)is,ISSTRIDE,&flg);
91: if (!flg) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_WRONG,"IS must be of type ISSTRIDE");
93: sub = (IS_Stride*)is->data;
94: if (first) *first = sub->first;
95: if (step) *step = sub->step;
96: return(0);
97: }
99: PetscErrorCode ISDestroy_Stride(IS is)
100: {
104: PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",NULL);
105: PetscFree(is->data);
106: return(0);
107: }
109: PetscErrorCode ISToGeneral_Stride(IS inis)
110: {
112: const PetscInt *idx;
113: PetscInt n;
116: ISGetLocalSize(inis,&n);
117: ISGetIndices(inis,&idx);
118: ISSetType(inis,ISGENERAL);
119: ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
120: return(0);
121: }
123: PetscErrorCode ISLocate_Stride(IS is,PetscInt key,PetscInt *location)
124: {
125: IS_Stride *sub = (IS_Stride*)is->data;
126: PetscInt rem, step;
129: *location = -1;
130: step = sub->step;
131: key -= sub->first;
132: rem = key / step;
133: if ((rem < is->map->n) && !(key % step)) {
134: *location = rem;
135: }
136: return(0);
137: }
139: /*
140: Returns a legitimate index memory even if
141: the stride index set is empty.
142: */
143: PetscErrorCode ISGetIndices_Stride(IS is,const PetscInt *idx[])
144: {
145: IS_Stride *sub = (IS_Stride*)is->data;
147: PetscInt i,**dx = (PetscInt**)idx;
150: PetscMalloc1(is->map->n,(PetscInt**)idx);
151: if (is->map->n) {
152: (*dx)[0] = sub->first;
153: for (i=1; i<is->map->n; i++) (*dx)[i] = (*dx)[i-1] + sub->step;
154: }
155: return(0);
156: }
158: PetscErrorCode ISRestoreIndices_Stride(IS in,const PetscInt *idx[])
159: {
163: PetscFree(*(void**)idx);
164: return(0);
165: }
167: PetscErrorCode ISView_Stride(IS is,PetscViewer viewer)
168: {
169: IS_Stride *sub = (IS_Stride*)is->data;
170: PetscInt i,n = is->map->n;
171: PetscMPIInt rank,size;
172: PetscBool iascii,ibinary;
173: PetscViewerFormat fmt;
174: PetscErrorCode ierr;
177: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
178: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&ibinary);
179: if (iascii) {
180: PetscBool matl, isperm;
182: MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
183: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
184: PetscViewerGetFormat(viewer,&fmt);
185: matl = (PetscBool)(fmt == PETSC_VIEWER_ASCII_MATLAB);
186: ISGetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_FALSE,&isperm);
187: if (isperm && !matl) {PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");}
188: if (size == 1) {
189: if (matl) {
190: const char* name;
192: PetscObjectGetName((PetscObject)is,&name);
193: PetscViewerASCIIPrintf(viewer,"%s = [%D : %D : %D];\n",name,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
194: } else {
195: PetscViewerASCIIPrintf(viewer,"Number of indices in (stride) set %D\n",n);
196: for (i=0; i<n; i++) {
197: PetscViewerASCIIPrintf(viewer,"%D %D\n",i,sub->first + i*sub->step);
198: }
199: }
200: PetscViewerFlush(viewer);
201: } else {
202: PetscViewerASCIIPushSynchronized(viewer);
203: if (matl) {
204: const char* name;
206: PetscObjectGetName((PetscObject)is,&name);
207: PetscViewerASCIISynchronizedPrintf(viewer,"%s_%d = [%D : %D : %D];\n",name,rank,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
208: } else {
209: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in (stride) set %D\n",rank,n);
210: for (i=0; i<n; i++) {
211: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,sub->first + i*sub->step);
212: }
213: }
214: PetscViewerFlush(viewer);
215: PetscViewerASCIIPopSynchronized(viewer);
216: }
217: } else if (ibinary) {
218: ISView_Binary(is,viewer);
219: }
220: return(0);
221: }
223: PetscErrorCode ISSort_Stride(IS is)
224: {
225: IS_Stride *sub = (IS_Stride*)is->data;
228: if (sub->step >= 0) return(0);
229: sub->first += (is->map->n - 1)*sub->step;
230: sub->step *= -1;
231: return(0);
232: }
234: PetscErrorCode ISSorted_Stride(IS is,PetscBool * flg)
235: {
236: IS_Stride *sub = (IS_Stride*)is->data;
239: if (sub->step >= 0) *flg = PETSC_TRUE;
240: else *flg = PETSC_FALSE;
241: return(0);
242: }
244: static PetscErrorCode ISUniqueLocal_Stride(IS is, PetscBool *flg)
245: {
246: IS_Stride *sub = (IS_Stride*)is->data;
249: if (!(is->map->n) || sub->step != 0) *flg = PETSC_TRUE;
250: else *flg = PETSC_FALSE;
251: return(0);
252: }
254: static PetscErrorCode ISPermutationLocal_Stride(IS is, PetscBool *flg)
255: {
256: IS_Stride *sub = (IS_Stride*)is->data;
259: if (!(is->map->n) || (PetscAbsInt(sub->step) == 1 && is->min == 0)) *flg = PETSC_TRUE;
260: else *flg = PETSC_FALSE;
261: return(0);
262: }
264: static PetscErrorCode ISIntervalLocal_Stride(IS is, PetscBool *flg)
265: {
266: IS_Stride *sub = (IS_Stride*)is->data;
269: if (!(is->map->n) || sub->step == 1) *flg = PETSC_TRUE;
270: else *flg = PETSC_FALSE;
271: return(0);
272: }
274: static PetscErrorCode ISOnComm_Stride(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
275: {
277: IS_Stride *sub = (IS_Stride*)is->data;
280: ISCreateStride(comm,is->map->n,sub->first,sub->step,newis);
281: return(0);
282: }
284: static PetscErrorCode ISSetBlockSize_Stride(IS is,PetscInt bs)
285: {
286: IS_Stride *sub = (IS_Stride*)is->data;
290: if (sub->step != 1 && bs != 1) SETERRQ2(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_SIZ,"ISSTRIDE has stride %D, cannot be blocked of size %D",sub->step,bs);
291: PetscLayoutSetBlockSize(is->map, bs);
292: return(0);
293: }
295: static PetscErrorCode ISContiguousLocal_Stride(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
296: {
297: IS_Stride *sub = (IS_Stride*)is->data;
300: if (sub->step == 1 && sub->first >= gstart && sub->first+is->map->n <= gend) {
301: *start = sub->first - gstart;
302: *contig = PETSC_TRUE;
303: } else {
304: *start = -1;
305: *contig = PETSC_FALSE;
306: }
307: return(0);
308: }
311: static struct _ISOps myops = { ISGetIndices_Stride,
312: ISRestoreIndices_Stride,
313: ISInvertPermutation_Stride,
314: ISSort_Stride,
315: ISSort_Stride,
316: ISSorted_Stride,
317: ISDuplicate_Stride,
318: ISDestroy_Stride,
319: ISView_Stride,
320: ISLoad_Default,
321: ISCopy_Stride,
322: ISToGeneral_Stride,
323: ISOnComm_Stride,
324: ISSetBlockSize_Stride,
325: ISContiguousLocal_Stride,
326: ISLocate_Stride,
327: ISSorted_Stride,
328: NULL,
329: ISUniqueLocal_Stride,
330: NULL,
331: ISPermutationLocal_Stride,
332: NULL,
333: ISIntervalLocal_Stride,
334: NULL};
337: /*@
338: ISStrideSetStride - Sets the stride information for a stride index set.
340: Collective on IS
342: Input Parameters:
343: + is - the index set
344: . n - the length of the locally owned portion of the index set
345: . first - the first element of the locally owned portion of the index set
346: - step - the change to the next index
348: Level: beginner
350: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather(), ISSTRIDE, ISCreateStride(), ISStrideGetInfo()
351: @*/
352: PetscErrorCode ISStrideSetStride(IS is,PetscInt n,PetscInt first,PetscInt step)
353: {
357: if (n < 0) SETERRQ1(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Negative length %D not valid", n);
358: ISClearInfoCache(is,PETSC_FALSE);
359: PetscUseMethod(is,"ISStrideSetStride_C",(IS,PetscInt,PetscInt,PetscInt),(is,n,first,step));
360: return(0);
361: }
363: PetscErrorCode ISStrideSetStride_Stride(IS is,PetscInt n,PetscInt first,PetscInt step)
364: {
366: PetscInt min,max;
367: IS_Stride *sub = (IS_Stride*)is->data;
368: PetscLayout map;
371: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is),n,is->map->N,is->map->bs,&map);
372: PetscLayoutDestroy(&is->map);
373: is->map = map;
375: sub->first = first;
376: sub->step = step;
377: if (step > 0) {min = first; max = first + step*(n-1);}
378: else {max = first; min = first + step*(n-1);}
380: is->min = n > 0 ? min : PETSC_MAX_INT;
381: is->max = n > 0 ? max : PETSC_MIN_INT;
382: is->data = (void*)sub;
383: return(0);
384: }
386: /*@
387: ISCreateStride - Creates a data structure for an index set
388: containing a list of evenly spaced integers.
390: Collective
392: Input Parameters:
393: + comm - the MPI communicator
394: . n - the length of the locally owned portion of the index set
395: . first - the first element of the locally owned portion of the index set
396: - step - the change to the next index
398: Output Parameter:
399: . is - the new index set
401: Notes:
402: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
403: conceptually the same as MPI_Group operations. The IS are the
404: distributed sets of indices and thus certain operations on them are collective.
406: Level: beginner
408: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather(), ISSTRIDE
409: @*/
410: PetscErrorCode ISCreateStride(MPI_Comm comm,PetscInt n,PetscInt first,PetscInt step,IS *is)
411: {
415: ISCreate(comm,is);
416: ISSetType(*is,ISSTRIDE);
417: ISStrideSetStride(*is,n,first,step);
418: return(0);
419: }
421: PETSC_EXTERN PetscErrorCode ISCreate_Stride(IS is)
422: {
424: IS_Stride *sub;
427: PetscNewLog(is,&sub);
428: is->data = (void *) sub;
429: PetscMemcpy(is->ops,&myops,sizeof(myops));
430: PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",ISStrideSetStride_Stride);
431: return(0);
432: }