Actual source code: stride.c
petsc-3.5.4 2015-05-23
2: /*
3: Index sets of evenly space integers, defined by a
4: start, stride and length.
5: */
6: #include <petsc-private/isimpl.h> /*I "petscis.h" I*/
7: #include <petscvec.h>
8: #include <petscviewer.h>
10: typedef struct {
11: PetscInt N,n,first,step;
12: } IS_Stride;
16: PetscErrorCode ISIdentity_Stride(IS is,PetscBool *ident)
17: {
18: IS_Stride *is_stride = (IS_Stride*)is->data;
21: is->isidentity = PETSC_FALSE;
22: *ident = PETSC_FALSE;
23: if (is_stride->first != 0) return(0);
24: if (is_stride->step != 1) return(0);
25: *ident = PETSC_TRUE;
26: is->isidentity = PETSC_TRUE;
27: return(0);
28: }
32: static PetscErrorCode ISCopy_Stride(IS is,IS isy)
33: {
34: IS_Stride *is_stride = (IS_Stride*)is->data,*isy_stride = (IS_Stride*)isy->data;
38: PetscMemcpy(isy_stride,is_stride,sizeof(IS_Stride));
39: return(0);
40: }
44: PetscErrorCode ISDuplicate_Stride(IS is,IS *newIS)
45: {
47: IS_Stride *sub = (IS_Stride*)is->data;
50: ISCreateStride(PetscObjectComm((PetscObject)is),sub->n,sub->first,sub->step,newIS);
51: return(0);
52: }
56: PetscErrorCode ISInvertPermutation_Stride(IS is,PetscInt nlocal,IS *perm)
57: {
58: IS_Stride *isstride = (IS_Stride*)is->data;
62: if (is->isidentity) {
63: ISCreateStride(PETSC_COMM_SELF,isstride->n,0,1,perm);
64: } else {
65: IS tmp;
66: const PetscInt *indices,n = isstride->n;
67: ISGetIndices(is,&indices);
68: ISCreateGeneral(PetscObjectComm((PetscObject)is),n,indices,PETSC_COPY_VALUES,&tmp);
69: ISSetPermutation(tmp);
70: ISRestoreIndices(is,&indices);
71: ISInvertPermutation(tmp,nlocal,perm);
72: ISDestroy(&tmp);
73: }
74: return(0);
75: }
79: /*@
80: ISStrideGetInfo - Returns the first index in a stride index set and
81: the stride width.
83: Not Collective
85: Input Parameter:
86: . is - the index set
88: Output Parameters:
89: . first - the first index
90: . step - the stride width
92: Level: intermediate
94: Notes:
95: Returns info on stride index set. This is a pseudo-public function that
96: should not be needed by most users.
98: Concepts: index sets^getting information
99: Concepts: IS^getting information
101: .seealso: ISCreateStride(), ISGetSize()
102: @*/
103: PetscErrorCode ISStrideGetInfo(IS is,PetscInt *first,PetscInt *step)
104: {
105: IS_Stride *sub;
106: PetscBool flg;
113: PetscObjectTypeCompare((PetscObject)is,ISSTRIDE,&flg);
114: if (!flg) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_WRONG,"IS must be of type ISSTRIDE");
116: sub = (IS_Stride*)is->data;
117: if (first) *first = sub->first;
118: if (step) *step = sub->step;
119: return(0);
120: }
124: PetscErrorCode ISDestroy_Stride(IS is)
125: {
129: PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",0);
130: PetscFree(is->data);
131: return(0);
132: }
136: PetscErrorCode ISToGeneral_Stride(IS inis)
137: {
139: const PetscInt *idx;
140: PetscInt n;
143: ISGetLocalSize(inis,&n);
144: ISGetIndices(inis,&idx);
145: ISSetType(inis,ISGENERAL);
146: ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
147: return(0);
148: }
151: /*
152: Returns a legitimate index memory even if
153: the stride index set is empty.
154: */
157: PetscErrorCode ISGetIndices_Stride(IS in,const PetscInt *idx[])
158: {
159: IS_Stride *sub = (IS_Stride*)in->data;
161: PetscInt i,**dx = (PetscInt**)idx;
164: PetscMalloc1(sub->n,idx);
165: if (sub->n) {
166: (*dx)[0] = sub->first;
167: for (i=1; i<sub->n; i++) (*dx)[i] = (*dx)[i-1] + sub->step;
168: }
169: return(0);
170: }
174: PetscErrorCode ISRestoreIndices_Stride(IS in,const PetscInt *idx[])
175: {
179: PetscFree(*(void**)idx);
180: return(0);
181: }
185: PetscErrorCode ISGetSize_Stride(IS is,PetscInt *size)
186: {
187: IS_Stride *sub = (IS_Stride*)is->data;
190: *size = sub->N;
191: return(0);
192: }
196: PetscErrorCode ISGetLocalSize_Stride(IS is,PetscInt *size)
197: {
198: IS_Stride *sub = (IS_Stride*)is->data;
201: *size = sub->n;
202: return(0);
203: }
207: PetscErrorCode ISView_Stride(IS is,PetscViewer viewer)
208: {
209: IS_Stride *sub = (IS_Stride*)is->data;
210: PetscInt i,n = sub->n;
211: PetscMPIInt rank,size;
212: PetscBool iascii;
216: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
217: if (iascii) {
218: MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
219: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
220: if (size == 1) {
221: if (is->isperm) {
222: PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");
223: }
224: PetscViewerASCIIPrintf(viewer,"Number of indices in (stride) set %D\n",n);
225: for (i=0; i<n; i++) {
226: PetscViewerASCIIPrintf(viewer,"%D %D\n",i,sub->first + i*sub->step);
227: }
228: PetscViewerFlush(viewer);
229: } else {
230: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
231: if (is->isperm) {
232: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutation\n",rank);
233: }
234: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in (stride) set %D\n",rank,n);
235: for (i=0; i<n; i++) {
236: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,sub->first + i*sub->step);
237: }
238: PetscViewerFlush(viewer);
239: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
240: }
241: }
242: return(0);
243: }
247: PetscErrorCode ISSort_Stride(IS is)
248: {
249: IS_Stride *sub = (IS_Stride*)is->data;
252: if (sub->step >= 0) return(0);
253: sub->first += (sub->n - 1)*sub->step;
254: sub->step *= -1;
255: return(0);
256: }
260: PetscErrorCode ISSorted_Stride(IS is,PetscBool * flg)
261: {
262: IS_Stride *sub = (IS_Stride*)is->data;
265: if (sub->step >= 0) *flg = PETSC_TRUE;
266: else *flg = PETSC_FALSE;
267: return(0);
268: }
272: static PetscErrorCode ISOnComm_Stride(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
273: {
275: IS_Stride *sub = (IS_Stride*)is->data;
278: ISCreateStride(comm,sub->n,sub->first,sub->step,newis);
279: return(0);
280: }
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: }
297: static PetscErrorCode ISContiguousLocal_Stride(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
298: {
299: IS_Stride *sub = (IS_Stride*)is->data;
302: if (sub->step == 1 && sub->first >= gstart && sub->first+sub->n <= gend) {
303: *start = sub->first - gstart;
304: *contig = PETSC_TRUE;
305: } else {
306: *start = -1;
307: *contig = PETSC_FALSE;
308: }
309: return(0);
310: }
313: static struct _ISOps myops = { ISGetSize_Stride,
314: ISGetLocalSize_Stride,
315: ISGetIndices_Stride,
316: ISRestoreIndices_Stride,
317: ISInvertPermutation_Stride,
318: ISSort_Stride,
319: ISSort_Stride,
320: ISSorted_Stride,
321: ISDuplicate_Stride,
322: ISDestroy_Stride,
323: ISView_Stride,
324: ISLoad_Default,
325: ISIdentity_Stride,
326: ISCopy_Stride,
327: ISToGeneral_Stride,
328: ISOnComm_Stride,
329: ISSetBlockSize_Stride,
330: ISContiguousLocal_Stride};
335: /*@
336: ISStrideSetStride - Sets the stride information for a stride index set.
338: Collective on IS
340: Input Parameters:
341: + is - the index set
342: . n - the length of the locally owned portion of the index set
343: . first - the first element of the locally owned portion of the index set
344: - step - the change to the next index
346: Level: beginner
348: Concepts: IS^stride
349: Concepts: index sets^stride
350: Concepts: stride^index set
352: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
353: @*/
354: PetscErrorCode ISStrideSetStride(IS is,PetscInt n,PetscInt first,PetscInt step)
355: {
359: if (n < 0) SETERRQ1(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Negative length %d not valid", n);
360: PetscUseMethod(is,"ISStrideSetStride_C",(IS,PetscInt,PetscInt,PetscInt),(is,n,first,step));
361: return(0);
362: }
366: PetscErrorCode ISStrideSetStride_Stride(IS is,PetscInt n,PetscInt first,PetscInt step)
367: {
369: PetscInt min,max;
370: IS_Stride *sub = (IS_Stride*)is->data;
373: sub->n = n;
374: MPI_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)is));
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 : -1;
381: is->max = n > 0 ? max : -1;
382: is->data = (void*)sub;
384: if ((!first && step == 1) || (first == max && step == -1 && !min)) is->isperm = PETSC_TRUE;
385: else is->isperm = PETSC_FALSE;
386: return(0);
387: }
391: /*@
392: ISCreateStride - Creates a data structure for an index set
393: containing a list of evenly spaced integers.
395: Collective on MPI_Comm
397: Input Parameters:
398: + comm - the MPI communicator
399: . n - the length of the locally owned portion of the index set
400: . first - the first element of the locally owned portion of the index set
401: - step - the change to the next index
403: Output Parameter:
404: . is - the new index set
406: Notes:
407: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
408: conceptually the same as MPI_Group operations. The IS are the
409: distributed sets of indices and thus certain operations on them are collective.
411: Level: beginner
413: Concepts: IS^stride
414: Concepts: index sets^stride
415: Concepts: stride^index set
417: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
418: @*/
419: PetscErrorCode ISCreateStride(MPI_Comm comm,PetscInt n,PetscInt first,PetscInt step,IS *is)
420: {
424: ISCreate(comm,is);
425: ISSetType(*is,ISSTRIDE);
426: ISStrideSetStride(*is,n,first,step);
427: return(0);
428: }
432: PETSC_EXTERN PetscErrorCode ISCreate_Stride(IS is)
433: {
435: IS_Stride *sub;
438: PetscMemcpy(is->ops,&myops,sizeof(myops));
439: PetscNewLog(is,&sub);
440: is->data = sub;
441: PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",ISStrideSetStride_Stride);
442: return(0);
443: }