Actual source code: stride.c
petsc-3.12.5 2020-03-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: PetscErrorCode ISIdentity_Stride(IS is,PetscBool *ident)
14: {
15: IS_Stride *is_stride = (IS_Stride*)is->data;
18: is->isidentity = PETSC_FALSE;
19: *ident = PETSC_FALSE;
20: if (is_stride->first != 0) return(0);
21: if (is_stride->step != 1) return(0);
22: *ident = PETSC_TRUE;
23: is->isidentity = PETSC_TRUE;
24: return(0);
25: }
27: static PetscErrorCode ISCopy_Stride(IS is,IS isy)
28: {
29: IS_Stride *is_stride = (IS_Stride*)is->data,*isy_stride = (IS_Stride*)isy->data;
33: PetscMemcpy(isy_stride,is_stride,sizeof(IS_Stride));
34: return(0);
35: }
37: PetscErrorCode ISDuplicate_Stride(IS is,IS *newIS)
38: {
40: IS_Stride *sub = (IS_Stride*)is->data;
43: ISCreateStride(PetscObjectComm((PetscObject)is),is->map->n,sub->first,sub->step,newIS);
44: return(0);
45: }
47: PetscErrorCode ISInvertPermutation_Stride(IS is,PetscInt nlocal,IS *perm)
48: {
52: if (is->isidentity) {
53: ISCreateStride(PETSC_COMM_SELF,is->map->n,0,1,perm);
54: } else {
55: IS tmp;
56: const PetscInt *indices,n = is->map->n;
57: ISGetIndices(is,&indices);
58: ISCreateGeneral(PetscObjectComm((PetscObject)is),n,indices,PETSC_COPY_VALUES,&tmp);
59: ISSetPermutation(tmp);
60: ISRestoreIndices(is,&indices);
61: ISInvertPermutation(tmp,nlocal,perm);
62: ISDestroy(&tmp);
63: }
64: return(0);
65: }
67: /*@
68: ISStrideGetInfo - Returns the first index in a stride index set and
69: the stride width.
71: Not Collective
73: Input Parameter:
74: . is - the index set
76: Output Parameters:
77: + first - the first index
78: - step - the stride width
80: Level: intermediate
82: Notes:
83: Returns info on stride index set. This is a pseudo-public function that
84: should not be needed by most users.
87: .seealso: ISCreateStride(), ISGetSize()
88: @*/
89: PetscErrorCode ISStrideGetInfo(IS is,PetscInt *first,PetscInt *step)
90: {
91: IS_Stride *sub;
92: PetscBool flg;
99: PetscObjectTypeCompare((PetscObject)is,ISSTRIDE,&flg);
100: if (!flg) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_WRONG,"IS must be of type ISSTRIDE");
102: sub = (IS_Stride*)is->data;
103: if (first) *first = sub->first;
104: if (step) *step = sub->step;
105: return(0);
106: }
108: PetscErrorCode ISDestroy_Stride(IS is)
109: {
113: PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",NULL);
114: PetscFree(is->data);
115: return(0);
116: }
118: PetscErrorCode ISToGeneral_Stride(IS inis)
119: {
121: const PetscInt *idx;
122: PetscInt n;
125: ISGetLocalSize(inis,&n);
126: ISGetIndices(inis,&idx);
127: ISSetType(inis,ISGENERAL);
128: ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
129: return(0);
130: }
132: PetscErrorCode ISLocate_Stride(IS is,PetscInt key,PetscInt *location)
133: {
134: IS_Stride *sub = (IS_Stride*)is->data;
135: PetscInt rem, step;
138: *location = -1;
139: step = sub->step;
140: key -= sub->first;
141: rem = key / step;
142: if ((rem < is->map->n) && !(key % step)) {
143: *location = rem;
144: }
145: return(0);
146: }
148: /*
149: Returns a legitimate index memory even if
150: the stride index set is empty.
151: */
152: PetscErrorCode ISGetIndices_Stride(IS is,const PetscInt *idx[])
153: {
154: IS_Stride *sub = (IS_Stride*)is->data;
156: PetscInt i,**dx = (PetscInt**)idx;
159: PetscMalloc1(is->map->n,(PetscInt**)idx);
160: if (is->map->n) {
161: (*dx)[0] = sub->first;
162: for (i=1; i<is->map->n; i++) (*dx)[i] = (*dx)[i-1] + sub->step;
163: }
164: return(0);
165: }
167: PetscErrorCode ISRestoreIndices_Stride(IS in,const PetscInt *idx[])
168: {
172: PetscFree(*(void**)idx);
173: return(0);
174: }
176: PetscErrorCode ISView_Stride(IS is,PetscViewer viewer)
177: {
178: IS_Stride *sub = (IS_Stride*)is->data;
179: PetscInt i,n = is->map->n;
180: PetscMPIInt rank,size;
181: PetscBool iascii;
182: PetscViewerFormat fmt;
183: PetscErrorCode ierr;
186: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
187: if (iascii) {
188: PetscBool matl;
190: MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
191: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
192: PetscViewerGetFormat(viewer,&fmt);
193: matl = (PetscBool)(fmt == PETSC_VIEWER_ASCII_MATLAB);
194: if (size == 1) {
195: if (matl) {
196: const char* name;
198: PetscObjectGetName((PetscObject)is,&name);
199: PetscViewerASCIIPrintf(viewer,"%s = [%D : %D : %D];\n",name,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
200: } else {
201: if (is->isperm) {
202: PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");
203: }
204: PetscViewerASCIIPrintf(viewer,"Number of indices in (stride) set %D\n",n);
205: for (i=0; i<n; i++) {
206: PetscViewerASCIIPrintf(viewer,"%D %D\n",i,sub->first + i*sub->step);
207: }
208: }
209: PetscViewerFlush(viewer);
210: } else {
211: PetscViewerASCIIPushSynchronized(viewer);
212: if (matl) {
213: const char* name;
215: PetscObjectGetName((PetscObject)is,&name);
216: PetscViewerASCIISynchronizedPrintf(viewer,"%s_%d = [%D : %D : %D];\n",name,rank,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
217: } else {
218: if (is->isperm) {
219: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutation\n",rank);
220: }
221: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in (stride) set %D\n",rank,n);
222: for (i=0; i<n; i++) {
223: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,sub->first + i*sub->step);
224: }
225: }
226: PetscViewerFlush(viewer);
227: PetscViewerASCIIPopSynchronized(viewer);
228: }
229: }
230: return(0);
231: }
233: PetscErrorCode ISSort_Stride(IS is)
234: {
235: IS_Stride *sub = (IS_Stride*)is->data;
238: if (sub->step >= 0) return(0);
239: sub->first += (is->map->n - 1)*sub->step;
240: sub->step *= -1;
241: return(0);
242: }
244: PetscErrorCode ISSorted_Stride(IS is,PetscBool * flg)
245: {
246: IS_Stride *sub = (IS_Stride*)is->data;
249: if (sub->step >= 0) *flg = PETSC_TRUE;
250: else *flg = PETSC_FALSE;
251: return(0);
252: }
254: static PetscErrorCode ISOnComm_Stride(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
255: {
257: IS_Stride *sub = (IS_Stride*)is->data;
260: ISCreateStride(comm,is->map->n,sub->first,sub->step,newis);
261: return(0);
262: }
264: static PetscErrorCode ISSetBlockSize_Stride(IS is,PetscInt bs)
265: {
266: IS_Stride *sub = (IS_Stride*)is->data;
270: 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);
271: PetscLayoutSetBlockSize(is->map, bs);
272: return(0);
273: }
275: static PetscErrorCode ISContiguousLocal_Stride(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
276: {
277: IS_Stride *sub = (IS_Stride*)is->data;
280: if (sub->step == 1 && sub->first >= gstart && sub->first+is->map->n <= gend) {
281: *start = sub->first - gstart;
282: *contig = PETSC_TRUE;
283: } else {
284: *start = -1;
285: *contig = PETSC_FALSE;
286: }
287: return(0);
288: }
291: static struct _ISOps myops = { ISGetIndices_Stride,
292: ISRestoreIndices_Stride,
293: ISInvertPermutation_Stride,
294: ISSort_Stride,
295: ISSort_Stride,
296: ISSorted_Stride,
297: ISDuplicate_Stride,
298: ISDestroy_Stride,
299: ISView_Stride,
300: ISLoad_Default,
301: ISIdentity_Stride,
302: ISCopy_Stride,
303: ISToGeneral_Stride,
304: ISOnComm_Stride,
305: ISSetBlockSize_Stride,
306: ISContiguousLocal_Stride,
307: ISLocate_Stride};
310: /*@
311: ISStrideSetStride - Sets the stride information for a stride index set.
313: Collective on IS
315: Input Parameters:
316: + is - the index set
317: . n - the length of the locally owned portion of the index set
318: . first - the first element of the locally owned portion of the index set
319: - step - the change to the next index
321: Level: beginner
323: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
324: @*/
325: PetscErrorCode ISStrideSetStride(IS is,PetscInt n,PetscInt first,PetscInt step)
326: {
330: if (n < 0) SETERRQ1(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Negative length %d not valid", n);
331: PetscUseMethod(is,"ISStrideSetStride_C",(IS,PetscInt,PetscInt,PetscInt),(is,n,first,step));
332: return(0);
333: }
335: PetscErrorCode ISStrideSetStride_Stride(IS is,PetscInt n,PetscInt first,PetscInt step)
336: {
338: PetscInt min,max;
339: IS_Stride *sub = (IS_Stride*)is->data;
340: PetscLayout map;
343: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is),n,is->map->N,is->map->bs,&map);
344: PetscLayoutDestroy(&is->map);
345: is->map = map;
347: sub->first = first;
348: sub->step = step;
349: if (step > 0) {min = first; max = first + step*(n-1);}
350: else {max = first; min = first + step*(n-1);}
352: is->min = n > 0 ? min : PETSC_MAX_INT;
353: is->max = n > 0 ? max : PETSC_MIN_INT;
354: is->data = (void*)sub;
356: if ((!first && step == 1) || (first == max && step == -1 && !min)) is->isperm = PETSC_TRUE;
357: else is->isperm = PETSC_FALSE;
358: is->isidentity = PETSC_FALSE;
359: return(0);
360: }
362: /*@
363: ISCreateStride - Creates a data structure for an index set
364: containing a list of evenly spaced integers.
366: Collective
368: Input Parameters:
369: + comm - the MPI communicator
370: . n - the length of the locally owned portion of the index set
371: . first - the first element of the locally owned portion of the index set
372: - step - the change to the next index
374: Output Parameter:
375: . is - the new index set
377: Notes:
378: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
379: conceptually the same as MPI_Group operations. The IS are the
380: distributed sets of indices and thus certain operations on them are collective.
382: Level: beginner
384: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
385: @*/
386: PetscErrorCode ISCreateStride(MPI_Comm comm,PetscInt n,PetscInt first,PetscInt step,IS *is)
387: {
391: ISCreate(comm,is);
392: ISSetType(*is,ISSTRIDE);
393: ISStrideSetStride(*is,n,first,step);
394: return(0);
395: }
397: PETSC_EXTERN PetscErrorCode ISCreate_Stride(IS is)
398: {
400: IS_Stride *sub;
403: PetscNewLog(is,&sub);
404: is->data = (void *) sub;
405: PetscMemcpy(is->ops,&myops,sizeof(myops));
406: PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",ISStrideSetStride_Stride);
407: return(0);
408: }