Actual source code: stride.c
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;
17: PetscMemcpy(isy_stride,is_stride,sizeof(IS_Stride));
18: return 0;
19: }
21: PetscErrorCode ISDuplicate_Stride(IS is,IS *newIS)
22: {
23: IS_Stride *sub = (IS_Stride*)is->data;
25: ISCreateStride(PetscObjectComm((PetscObject)is),is->map->n,sub->first,sub->step,newIS);
26: return 0;
27: }
29: PetscErrorCode ISInvertPermutation_Stride(IS is,PetscInt nlocal,IS *perm)
30: {
31: PetscBool isident;
33: ISGetInfo(is,IS_IDENTITY,IS_GLOBAL,PETSC_TRUE,&isident);
34: if (isident) {
35: PetscInt rStart, rEnd;
37: PetscLayoutGetRange(is->map, &rStart, &rEnd);
38: ISCreateStride(PETSC_COMM_SELF,PetscMax(rEnd - rStart, 0),rStart,1,perm);
39: } else {
40: IS tmp;
41: const PetscInt *indices,n = is->map->n;
43: ISGetIndices(is,&indices);
44: ISCreateGeneral(PetscObjectComm((PetscObject)is),n,indices,PETSC_COPY_VALUES,&tmp);
45: ISSetPermutation(tmp);
46: ISRestoreIndices(is,&indices);
47: ISInvertPermutation(tmp,nlocal,perm);
48: ISDestroy(&tmp);
49: }
50: return 0;
51: }
53: /*@
54: ISStrideGetInfo - Returns the first index in a stride index set and the stride width.
56: Not Collective
58: Input Parameter:
59: . is - the index set
61: Output Parameters:
62: + first - the first index
63: - step - the stride width
65: Level: intermediate
67: Notes:
68: Returns info on stride index set. This is a pseudo-public function that
69: should not be needed by most users.
71: .seealso: ISCreateStride(), ISGetSize(), ISSTRIDE
72: @*/
73: PetscErrorCode ISStrideGetInfo(IS is,PetscInt *first,PetscInt *step)
74: {
75: IS_Stride *sub;
76: PetscBool flg;
81: PetscObjectTypeCompare((PetscObject)is,ISSTRIDE,&flg);
84: sub = (IS_Stride*)is->data;
85: if (first) *first = sub->first;
86: if (step) *step = sub->step;
87: return 0;
88: }
90: PetscErrorCode ISDestroy_Stride(IS is)
91: {
92: PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",NULL);
93: PetscFree(is->data);
94: return 0;
95: }
97: PetscErrorCode ISToGeneral_Stride(IS inis)
98: {
99: const PetscInt *idx;
100: PetscInt n;
102: ISGetLocalSize(inis,&n);
103: ISGetIndices(inis,&idx);
104: ISSetType(inis,ISGENERAL);
105: ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
106: return 0;
107: }
109: PetscErrorCode ISLocate_Stride(IS is,PetscInt key,PetscInt *location)
110: {
111: IS_Stride *sub = (IS_Stride*)is->data;
112: PetscInt rem, step;
114: *location = -1;
115: step = sub->step;
116: key -= sub->first;
117: rem = key / step;
118: if ((rem < is->map->n) && !(key % step)) {
119: *location = rem;
120: }
121: return 0;
122: }
124: /*
125: Returns a legitimate index memory even if
126: the stride index set is empty.
127: */
128: PetscErrorCode ISGetIndices_Stride(IS is,const PetscInt *idx[])
129: {
130: IS_Stride *sub = (IS_Stride*)is->data;
131: PetscInt i,**dx = (PetscInt**)idx;
133: PetscMalloc1(is->map->n,(PetscInt**)idx);
134: if (is->map->n) {
135: (*dx)[0] = sub->first;
136: for (i=1; i<is->map->n; i++) (*dx)[i] = (*dx)[i-1] + sub->step;
137: }
138: return 0;
139: }
141: PetscErrorCode ISRestoreIndices_Stride(IS in,const PetscInt *idx[])
142: {
143: PetscFree(*(void**)idx);
144: return 0;
145: }
147: PetscErrorCode ISView_Stride(IS is,PetscViewer viewer)
148: {
149: IS_Stride *sub = (IS_Stride*)is->data;
150: PetscInt i,n = is->map->n;
151: PetscMPIInt rank,size;
152: PetscBool iascii,ibinary;
153: PetscViewerFormat fmt;
155: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
156: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&ibinary);
157: if (iascii) {
158: PetscBool matl, isperm;
160: MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
161: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
162: PetscViewerGetFormat(viewer,&fmt);
163: matl = (PetscBool)(fmt == PETSC_VIEWER_ASCII_MATLAB);
164: ISGetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_FALSE,&isperm);
165: if (isperm && !matl) PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");
166: if (size == 1) {
167: if (matl) {
168: const char* name;
170: PetscObjectGetName((PetscObject)is,&name);
171: PetscViewerASCIIPrintf(viewer,"%s = [%" PetscInt_FMT " : %" PetscInt_FMT " : %" PetscInt_FMT "];\n",name,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
172: } else {
173: PetscViewerASCIIPrintf(viewer,"Number of indices in (stride) set %" PetscInt_FMT "\n",n);
174: for (i=0; i<n; i++) {
175: PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT "\n",i,sub->first + i*sub->step);
176: }
177: }
178: PetscViewerFlush(viewer);
179: } else {
180: PetscViewerASCIIPushSynchronized(viewer);
181: if (matl) {
182: const char* name;
184: PetscObjectGetName((PetscObject)is,&name);
185: PetscViewerASCIISynchronizedPrintf(viewer,"%s_%d = [%" PetscInt_FMT " : %" PetscInt_FMT " : %" PetscInt_FMT "];\n",name,rank,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
186: } else {
187: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in (stride) set %" PetscInt_FMT "\n",rank,n);
188: for (i=0; i<n; i++) {
189: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %" PetscInt_FMT " %" PetscInt_FMT "\n",rank,i,sub->first + i*sub->step);
190: }
191: }
192: PetscViewerFlush(viewer);
193: PetscViewerASCIIPopSynchronized(viewer);
194: }
195: } else if (ibinary) {
196: ISView_Binary(is,viewer);
197: }
198: return 0;
199: }
201: PetscErrorCode ISSort_Stride(IS is)
202: {
203: IS_Stride *sub = (IS_Stride*)is->data;
205: if (sub->step >= 0) return 0;
206: sub->first += (is->map->n - 1)*sub->step;
207: sub->step *= -1;
208: return 0;
209: }
211: PetscErrorCode ISSorted_Stride(IS is,PetscBool * flg)
212: {
213: IS_Stride *sub = (IS_Stride*)is->data;
215: if (sub->step >= 0) *flg = PETSC_TRUE;
216: else *flg = PETSC_FALSE;
217: return 0;
218: }
220: static PetscErrorCode ISUniqueLocal_Stride(IS is, PetscBool *flg)
221: {
222: IS_Stride *sub = (IS_Stride*)is->data;
224: if (!(is->map->n) || sub->step != 0) *flg = PETSC_TRUE;
225: else *flg = PETSC_FALSE;
226: return 0;
227: }
229: static PetscErrorCode ISPermutationLocal_Stride(IS is, PetscBool *flg)
230: {
231: IS_Stride *sub = (IS_Stride*)is->data;
233: if (!(is->map->n) || (PetscAbsInt(sub->step) == 1 && is->min == 0)) *flg = PETSC_TRUE;
234: else *flg = PETSC_FALSE;
235: return 0;
236: }
238: static PetscErrorCode ISIntervalLocal_Stride(IS is, PetscBool *flg)
239: {
240: IS_Stride *sub = (IS_Stride*)is->data;
242: if (!(is->map->n) || sub->step == 1) *flg = PETSC_TRUE;
243: else *flg = PETSC_FALSE;
244: return 0;
245: }
247: static PetscErrorCode ISOnComm_Stride(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
248: {
249: IS_Stride *sub = (IS_Stride*)is->data;
251: ISCreateStride(comm,is->map->n,sub->first,sub->step,newis);
252: return 0;
253: }
255: static PetscErrorCode ISSetBlockSize_Stride(IS is,PetscInt bs)
256: {
257: IS_Stride *sub = (IS_Stride*)is->data;
260: PetscLayoutSetBlockSize(is->map, bs);
261: return 0;
262: }
264: static PetscErrorCode ISContiguousLocal_Stride(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
265: {
266: IS_Stride *sub = (IS_Stride*)is->data;
268: if (sub->step == 1 && sub->first >= gstart && sub->first+is->map->n <= gend) {
269: *start = sub->first - gstart;
270: *contig = PETSC_TRUE;
271: } else {
272: *start = -1;
273: *contig = PETSC_FALSE;
274: }
275: return 0;
276: }
278: static struct _ISOps myops = {
279: PetscDesignatedInitializer(getindices,ISGetIndices_Stride),
280: PetscDesignatedInitializer(restoreindices,ISRestoreIndices_Stride),
281: PetscDesignatedInitializer(invertpermutation,ISInvertPermutation_Stride),
282: PetscDesignatedInitializer(sort,ISSort_Stride),
283: PetscDesignatedInitializer(sortremovedups,ISSort_Stride),
284: PetscDesignatedInitializer(sorted,ISSorted_Stride),
285: PetscDesignatedInitializer(duplicate,ISDuplicate_Stride),
286: PetscDesignatedInitializer(destroy,ISDestroy_Stride),
287: PetscDesignatedInitializer(view,ISView_Stride),
288: PetscDesignatedInitializer(load,ISLoad_Default),
289: PetscDesignatedInitializer(copy,ISCopy_Stride),
290: PetscDesignatedInitializer(togeneral,ISToGeneral_Stride),
291: PetscDesignatedInitializer(oncomm,ISOnComm_Stride),
292: PetscDesignatedInitializer(setblocksize,ISSetBlockSize_Stride),
293: PetscDesignatedInitializer(contiguous,ISContiguousLocal_Stride),
294: PetscDesignatedInitializer(locate,ISLocate_Stride),
295: PetscDesignatedInitializer(sortedlocal,ISSorted_Stride),
296: NULL,
297: PetscDesignatedInitializer(uniquelocal,ISUniqueLocal_Stride),
298: NULL,
299: PetscDesignatedInitializer(permlocal,ISPermutationLocal_Stride),
300: NULL,
301: PetscDesignatedInitializer(intervallocal,ISIntervalLocal_Stride),
302: NULL
303: };
305: /*@
306: ISStrideSetStride - Sets the stride information for a stride index set.
308: Collective on IS
310: Input Parameters:
311: + is - the index set
312: . n - the length of the locally owned portion of the index set
313: . first - the first element of the locally owned portion of the index set
314: - step - the change to the next index
316: Level: beginner
318: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather(), ISSTRIDE, ISCreateStride(), ISStrideGetInfo()
319: @*/
320: PetscErrorCode ISStrideSetStride(IS is,PetscInt n,PetscInt first,PetscInt step)
321: {
323: ISClearInfoCache(is,PETSC_FALSE);
324: PetscUseMethod(is,"ISStrideSetStride_C",(IS,PetscInt,PetscInt,PetscInt),(is,n,first,step));
325: return 0;
326: }
328: PetscErrorCode ISStrideSetStride_Stride(IS is,PetscInt n,PetscInt first,PetscInt step)
329: {
330: PetscInt min,max;
331: IS_Stride *sub = (IS_Stride*)is->data;
332: PetscLayout map;
334: PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is),n,is->map->N,is->map->bs,&map);
335: PetscLayoutDestroy(&is->map);
336: is->map = map;
338: sub->first = first;
339: sub->step = step;
340: if (step > 0) {min = first; max = first + step*(n-1);}
341: else {max = first; min = first + step*(n-1);}
343: is->min = n > 0 ? min : PETSC_MAX_INT;
344: is->max = n > 0 ? max : PETSC_MIN_INT;
345: is->data = (void*)sub;
346: return 0;
347: }
349: /*@
350: ISCreateStride - Creates a data structure for an index set
351: containing a list of evenly spaced integers.
353: Collective
355: Input Parameters:
356: + comm - the MPI communicator
357: . n - the length of the locally owned portion of the index set
358: . first - the first element of the locally owned portion of the index set
359: - step - the change to the next index
361: Output Parameter:
362: . is - the new index set
364: Notes:
365: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
366: conceptually the same as MPI_Group operations. The IS are the
367: distributed sets of indices and thus certain operations on them are collective.
369: Level: beginner
371: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather(), ISSTRIDE
372: @*/
373: PetscErrorCode ISCreateStride(MPI_Comm comm,PetscInt n,PetscInt first,PetscInt step,IS *is)
374: {
375: ISCreate(comm,is);
376: ISSetType(*is,ISSTRIDE);
377: ISStrideSetStride(*is,n,first,step);
378: return 0;
379: }
381: PETSC_EXTERN PetscErrorCode ISCreate_Stride(IS is)
382: {
383: IS_Stride *sub;
385: PetscNewLog(is,&sub);
386: is->data = (void *) sub;
387: PetscMemcpy(is->ops,&myops,sizeof(myops));
388: PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",ISStrideSetStride_Stride);
389: return 0;
390: }