Actual source code: stride.c
petsc-3.13.6 2020-09-29
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
61: the stride width.
63: Not Collective
65: Input Parameter:
66: . is - the index set
68: Output Parameters:
69: + first - the first index
70: - step - the stride width
72: Level: intermediate
74: Notes:
75: Returns info on stride index set. This is a pseudo-public function that
76: should not be needed by most users.
79: .seealso: ISCreateStride(), ISGetSize()
80: @*/
81: PetscErrorCode ISStrideGetInfo(IS is,PetscInt *first,PetscInt *step)
82: {
83: IS_Stride *sub;
84: PetscBool flg;
91: PetscObjectTypeCompare((PetscObject)is,ISSTRIDE,&flg);
92: if (!flg) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_WRONG,"IS must be of type ISSTRIDE");
94: sub = (IS_Stride*)is->data;
95: if (first) *first = sub->first;
96: if (step) *step = sub->step;
97: return(0);
98: }
100: PetscErrorCode ISDestroy_Stride(IS is)
101: {
105: PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",NULL);
106: PetscFree(is->data);
107: return(0);
108: }
110: PetscErrorCode ISToGeneral_Stride(IS inis)
111: {
113: const PetscInt *idx;
114: PetscInt n;
117: ISGetLocalSize(inis,&n);
118: ISGetIndices(inis,&idx);
119: ISSetType(inis,ISGENERAL);
120: ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
121: return(0);
122: }
124: PetscErrorCode ISLocate_Stride(IS is,PetscInt key,PetscInt *location)
125: {
126: IS_Stride *sub = (IS_Stride*)is->data;
127: PetscInt rem, step;
130: *location = -1;
131: step = sub->step;
132: key -= sub->first;
133: rem = key / step;
134: if ((rem < is->map->n) && !(key % step)) {
135: *location = rem;
136: }
137: return(0);
138: }
140: /*
141: Returns a legitimate index memory even if
142: the stride index set is empty.
143: */
144: PetscErrorCode ISGetIndices_Stride(IS is,const PetscInt *idx[])
145: {
146: IS_Stride *sub = (IS_Stride*)is->data;
148: PetscInt i,**dx = (PetscInt**)idx;
151: PetscMalloc1(is->map->n,(PetscInt**)idx);
152: if (is->map->n) {
153: (*dx)[0] = sub->first;
154: for (i=1; i<is->map->n; i++) (*dx)[i] = (*dx)[i-1] + sub->step;
155: }
156: return(0);
157: }
159: PetscErrorCode ISRestoreIndices_Stride(IS in,const PetscInt *idx[])
160: {
164: PetscFree(*(void**)idx);
165: return(0);
166: }
168: PetscErrorCode ISView_Stride(IS is,PetscViewer viewer)
169: {
170: IS_Stride *sub = (IS_Stride*)is->data;
171: PetscInt i,n = is->map->n;
172: PetscMPIInt rank,size;
173: PetscBool iascii,ibinary;
174: PetscViewerFormat fmt;
175: PetscErrorCode ierr;
178: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
179: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&ibinary);
180: if (iascii) {
181: PetscBool matl, isperm;
183: MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
184: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
185: PetscViewerGetFormat(viewer,&fmt);
186: matl = (PetscBool)(fmt == PETSC_VIEWER_ASCII_MATLAB);
187: ISGetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_FALSE,&isperm);
188: if (isperm && !matl) {PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");}
189: if (size == 1) {
190: if (matl) {
191: const char* name;
193: PetscObjectGetName((PetscObject)is,&name);
194: PetscViewerASCIIPrintf(viewer,"%s = [%D : %D : %D];\n",name,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
195: } else {
196: PetscViewerASCIIPrintf(viewer,"Number of indices in (stride) set %D\n",n);
197: for (i=0; i<n; i++) {
198: PetscViewerASCIIPrintf(viewer,"%D %D\n",i,sub->first + i*sub->step);
199: }
200: }
201: PetscViewerFlush(viewer);
202: } else {
203: PetscViewerASCIIPushSynchronized(viewer);
204: if (matl) {
205: const char* name;
207: PetscObjectGetName((PetscObject)is,&name);
208: PetscViewerASCIISynchronizedPrintf(viewer,"%s_%d = [%D : %D : %D];\n",name,rank,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
209: } else {
210: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in (stride) set %D\n",rank,n);
211: for (i=0; i<n; i++) {
212: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,sub->first + i*sub->step);
213: }
214: }
215: PetscViewerFlush(viewer);
216: PetscViewerASCIIPopSynchronized(viewer);
217: }
218: } else if (ibinary) {
219: ISView_Binary(is,viewer);
220: }
221: return(0);
222: }
224: PetscErrorCode ISSort_Stride(IS is)
225: {
226: IS_Stride *sub = (IS_Stride*)is->data;
229: if (sub->step >= 0) return(0);
230: sub->first += (is->map->n - 1)*sub->step;
231: sub->step *= -1;
232: return(0);
233: }
235: PetscErrorCode ISSorted_Stride(IS is,PetscBool * flg)
236: {
237: IS_Stride *sub = (IS_Stride*)is->data;
240: if (sub->step >= 0) *flg = PETSC_TRUE;
241: else *flg = PETSC_FALSE;
242: return(0);
243: }
245: static PetscErrorCode ISUniqueLocal_Stride(IS is, PetscBool *flg)
246: {
247: IS_Stride *sub = (IS_Stride*)is->data;
250: if (!(is->map->n) || sub->step != 0) *flg = PETSC_TRUE;
251: else *flg = PETSC_FALSE;
252: return(0);
253: }
255: static PetscErrorCode ISPermutationLocal_Stride(IS is, PetscBool *flg)
256: {
257: IS_Stride *sub = (IS_Stride*)is->data;
260: if (!(is->map->n) || (PetscAbsInt(sub->step) == 1 && is->min == 0)) *flg = PETSC_TRUE;
261: else *flg = PETSC_FALSE;
262: return(0);
263: }
265: static PetscErrorCode ISIntervalLocal_Stride(IS is, PetscBool *flg)
266: {
267: IS_Stride *sub = (IS_Stride*)is->data;
270: if (!(is->map->n) || sub->step == 1) *flg = PETSC_TRUE;
271: else *flg = PETSC_FALSE;
272: return(0);
273: }
275: static PetscErrorCode ISOnComm_Stride(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
276: {
278: IS_Stride *sub = (IS_Stride*)is->data;
281: ISCreateStride(comm,is->map->n,sub->first,sub->step,newis);
282: return(0);
283: }
285: static PetscErrorCode ISSetBlockSize_Stride(IS is,PetscInt bs)
286: {
287: IS_Stride *sub = (IS_Stride*)is->data;
291: 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);
292: PetscLayoutSetBlockSize(is->map, bs);
293: return(0);
294: }
296: static PetscErrorCode ISContiguousLocal_Stride(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
297: {
298: IS_Stride *sub = (IS_Stride*)is->data;
301: if (sub->step == 1 && sub->first >= gstart && sub->first+is->map->n <= gend) {
302: *start = sub->first - gstart;
303: *contig = PETSC_TRUE;
304: } else {
305: *start = -1;
306: *contig = PETSC_FALSE;
307: }
308: return(0);
309: }
312: static struct _ISOps myops = { ISGetIndices_Stride,
313: ISRestoreIndices_Stride,
314: ISInvertPermutation_Stride,
315: ISSort_Stride,
316: ISSort_Stride,
317: ISSorted_Stride,
318: ISDuplicate_Stride,
319: ISDestroy_Stride,
320: ISView_Stride,
321: ISLoad_Default,
322: ISCopy_Stride,
323: ISToGeneral_Stride,
324: ISOnComm_Stride,
325: ISSetBlockSize_Stride,
326: ISContiguousLocal_Stride,
327: ISLocate_Stride,
328: ISSorted_Stride,
329: NULL,
330: ISUniqueLocal_Stride,
331: NULL,
332: ISPermutationLocal_Stride,
333: NULL,
334: ISIntervalLocal_Stride,
335: NULL};
338: /*@
339: ISStrideSetStride - Sets the stride information for a stride index set.
341: Collective on IS
343: Input Parameters:
344: + is - the index set
345: . n - the length of the locally owned portion of the index set
346: . first - the first element of the locally owned portion of the index set
347: - step - the change to the next index
349: Level: beginner
351: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
352: @*/
353: PetscErrorCode ISStrideSetStride(IS is,PetscInt n,PetscInt first,PetscInt step)
354: {
358: if (n < 0) SETERRQ1(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Negative length %D not valid", n);
359: ISClearInfoCache(is,PETSC_FALSE);
360: PetscUseMethod(is,"ISStrideSetStride_C",(IS,PetscInt,PetscInt,PetscInt),(is,n,first,step));
361: return(0);
362: }
364: PetscErrorCode ISStrideSetStride_Stride(IS is,PetscInt n,PetscInt first,PetscInt step)
365: {
367: PetscInt min,max;
368: IS_Stride *sub = (IS_Stride*)is->data;
369: PetscLayout map;
372: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is),n,is->map->N,is->map->bs,&map);
373: PetscLayoutDestroy(&is->map);
374: is->map = map;
376: sub->first = first;
377: sub->step = step;
378: if (step > 0) {min = first; max = first + step*(n-1);}
379: else {max = first; min = first + step*(n-1);}
381: is->min = n > 0 ? min : PETSC_MAX_INT;
382: is->max = n > 0 ? max : PETSC_MIN_INT;
383: is->data = (void*)sub;
384: return(0);
385: }
387: /*@
388: ISCreateStride - Creates a data structure for an index set
389: containing a list of evenly spaced integers.
391: Collective
393: Input Parameters:
394: + comm - the MPI communicator
395: . n - the length of the locally owned portion of the index set
396: . first - the first element of the locally owned portion of the index set
397: - step - the change to the next index
399: Output Parameter:
400: . is - the new index set
402: Notes:
403: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
404: conceptually the same as MPI_Group operations. The IS are the
405: distributed sets of indices and thus certain operations on them are collective.
407: Level: beginner
409: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
410: @*/
411: PetscErrorCode ISCreateStride(MPI_Comm comm,PetscInt n,PetscInt first,PetscInt step,IS *is)
412: {
416: ISCreate(comm,is);
417: ISSetType(*is,ISSTRIDE);
418: ISStrideSetStride(*is,n,first,step);
419: return(0);
420: }
422: PETSC_EXTERN PetscErrorCode ISCreate_Stride(IS is)
423: {
425: IS_Stride *sub;
428: PetscNewLog(is,&sub);
429: is->data = (void *) sub;
430: PetscMemcpy(is->ops,&myops,sizeof(myops));
431: PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",ISStrideSetStride_Stride);
432: return(0);
433: }