Actual source code: stride.c
petsc-3.3-p7 2013-05-11
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>
9: typedef struct {
10: PetscInt N,n,first,step;
11: } IS_Stride;
15: PetscErrorCode ISIdentity_Stride(IS is,PetscBool *ident)
16: {
17: IS_Stride *is_stride = (IS_Stride*)is->data;
20: is->isidentity = PETSC_FALSE;
21: *ident = PETSC_FALSE;
22: if (is_stride->first != 0) return(0);
23: if (is_stride->step != 1) return(0);
24: *ident = PETSC_TRUE;
25: is->isidentity = PETSC_TRUE;
26: return(0);
27: }
31: static PetscErrorCode ISCopy_Stride(IS is,IS isy)
32: {
33: IS_Stride *is_stride = (IS_Stride*)is->data,*isy_stride = (IS_Stride*)isy->data;
37: PetscMemcpy(isy_stride,is_stride,sizeof(IS_Stride));
38: return(0);
39: }
43: PetscErrorCode ISDuplicate_Stride(IS is,IS *newIS)
44: {
46: IS_Stride *sub = (IS_Stride*)is->data;
49: ISCreateStride(((PetscObject)is)->comm,sub->n,sub->first,sub->step,newIS);
50: return(0);
51: }
55: PetscErrorCode ISInvertPermutation_Stride(IS is,PetscInt nlocal,IS *perm)
56: {
57: IS_Stride *isstride = (IS_Stride*)is->data;
61: if (is->isidentity) {
62: ISCreateStride(PETSC_COMM_SELF,isstride->n,0,1,perm);
63: } else {
64: IS tmp;
65: const PetscInt *indices,n = isstride->n;
66: ISGetIndices(is,&indices);
67: ISCreateGeneral(((PetscObject)is)->comm,n,indices,PETSC_COPY_VALUES,&tmp);
68: ISSetPermutation(tmp);
69: ISRestoreIndices(is,&indices);
70: ISInvertPermutation(tmp,nlocal,perm);
71: ISDestroy(&tmp);
72: }
73: return(0);
74: }
75:
78: /*@
79: ISStrideGetInfo - Returns the first index in a stride index set and
80: the stride width.
82: Not Collective
84: Input Parameter:
85: . is - the index set
87: Output Parameters:
88: . first - the first index
89: . step - the stride width
91: Level: intermediate
93: Notes:
94: Returns info on stride index set. This is a pseudo-public function that
95: should not be needed by most users.
97: Concepts: index sets^getting information
98: Concepts: IS^getting information
100: .seealso: ISCreateStride(), ISGetSize()
101: @*/
102: PetscErrorCode ISStrideGetInfo(IS is,PetscInt *first,PetscInt *step)
103: {
104: IS_Stride *sub;
111: sub = (IS_Stride*)is->data;
112: if (first) *first = sub->first;
113: if (step) *step = sub->step;
114: return(0);
115: }
119: PetscErrorCode ISDestroy_Stride(IS is)
120: {
124: PetscObjectComposeFunctionDynamic((PetscObject)is,"ISStrideSetStride_C","",0);
125: PetscFree(is->data);
126: return(0);
127: }
131: PetscErrorCode ISToGeneral_Stride(IS inis)
132: {
134: const PetscInt *idx;
135: PetscInt n;
138: ISGetLocalSize(inis,&n);
139: ISGetIndices(inis,&idx);
140: ISSetType(inis,ISGENERAL);
141: ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
142: return(0);
143: }
146: /*
147: Returns a legitimate index memory even if
148: the stride index set is empty.
149: */
152: PetscErrorCode ISGetIndices_Stride(IS in,const PetscInt *idx[])
153: {
154: IS_Stride *sub = (IS_Stride*)in->data;
156: PetscInt i,**dx = (PetscInt**)idx;
159: PetscMalloc(sub->n*sizeof(PetscInt),idx);
160: if (sub->n) {
161: (*dx)[0] = sub->first;
162: for (i=1; i<sub->n; i++) (*dx)[i] = (*dx)[i-1] + sub->step;
163: }
164: return(0);
165: }
169: PetscErrorCode ISRestoreIndices_Stride(IS in,const PetscInt *idx[])
170: {
174: PetscFree(*(void**)idx);
175: return(0);
176: }
180: PetscErrorCode ISGetSize_Stride(IS is,PetscInt *size)
181: {
182: IS_Stride *sub = (IS_Stride *)is->data;
185: *size = sub->N;
186: return(0);
187: }
191: PetscErrorCode ISGetLocalSize_Stride(IS is,PetscInt *size)
192: {
193: IS_Stride *sub = (IS_Stride *)is->data;
196: *size = sub->n;
197: return(0);
198: }
202: PetscErrorCode ISView_Stride(IS is,PetscViewer viewer)
203: {
204: IS_Stride *sub = (IS_Stride *)is->data;
205: PetscInt i,n = sub->n;
206: PetscMPIInt rank,size;
207: PetscBool iascii;
211: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
212: if (iascii) {
213: MPI_Comm_rank(((PetscObject)is)->comm,&rank);
214: MPI_Comm_size(((PetscObject)is)->comm,&size);
215: if (size == 1) {
216: if (is->isperm) {
217: PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");
218: }
219: PetscViewerASCIIPrintf(viewer,"Number of indices in (stride) set %D\n",n);
220: for (i=0; i<n; i++) {
221: PetscViewerASCIIPrintf(viewer,"%D %D\n",i,sub->first + i*sub->step);
222: }
223: PetscViewerFlush(viewer);
224: } else {
225: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
226: if (is->isperm) {
227: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutation\n",rank);
228: }
229: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in (stride) set %D\n",rank,n);
230: for (i=0; i<n; i++) {
231: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,sub->first + i*sub->step);
232: }
233: PetscViewerFlush(viewer);
234: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
235: }
236: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for this object",((PetscObject)viewer)->type_name);
237: return(0);
238: }
239:
242: PetscErrorCode ISSort_Stride(IS is)
243: {
244: IS_Stride *sub = (IS_Stride*)is->data;
247: if (sub->step >= 0) return(0);
248: sub->first += (sub->n - 1)*sub->step;
249: sub->step *= -1;
250: return(0);
251: }
255: PetscErrorCode ISSorted_Stride(IS is,PetscBool * flg)
256: {
257: IS_Stride *sub = (IS_Stride*)is->data;
260: if (sub->step >= 0) *flg = PETSC_TRUE;
261: else *flg = PETSC_FALSE;
262: return(0);
263: }
267: static PetscErrorCode ISOnComm_Stride(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
268: {
270: IS_Stride *sub = (IS_Stride*)is->data;
273: ISCreateStride(comm,sub->n,sub->first,sub->step,newis);
274: return(0);
275: }
279: static PetscErrorCode ISSetBlockSize_Stride(IS is,PetscInt bs)
280: {
281: IS_Stride *sub = (IS_Stride*)is->data;
284: if (sub->step != 1 && bs != 1) SETERRQ2(((PetscObject)is)->comm,PETSC_ERR_ARG_SIZ,"ISSTRIDE has stride %D, cannot be blocked of size %D",sub->step,bs);
285: is->bs = bs;
286: return(0);
287: }
291: static PetscErrorCode ISContiguousLocal_Stride(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
292: {
293: IS_Stride *sub = (IS_Stride*)is->data;
296: if (sub->step == 1 && sub->first >= gstart && sub->first+sub->n <= gend) {
297: *start = sub->first - gstart;
298: *contig = PETSC_TRUE;
299: } else {
300: *start = -1;
301: *contig = PETSC_FALSE;
302: }
303: return(0);
304: }
307: static struct _ISOps myops = { ISGetSize_Stride,
308: ISGetLocalSize_Stride,
309: ISGetIndices_Stride,
310: ISRestoreIndices_Stride,
311: ISInvertPermutation_Stride,
312: ISSort_Stride,
313: ISSorted_Stride,
314: ISDuplicate_Stride,
315: ISDestroy_Stride,
316: ISView_Stride,
317: ISIdentity_Stride,
318: ISCopy_Stride,
319: ISToGeneral_Stride,
320: ISOnComm_Stride,
321: ISSetBlockSize_Stride,
322: ISContiguousLocal_Stride
323: };
328: /*@
329: ISStrideSetStride - Sets the stride information for a stride index set.
331: Collective on IS
333: Input Parameters:
334: + is - the index set
335: . n - the length of the locally owned portion of the index set
336: . first - the first element of the locally owned portion of the index set
337: - step - the change to the next index
339: Level: beginner
341: Concepts: IS^stride
342: Concepts: index sets^stride
343: Concepts: stride^index set
345: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
346: @*/
347: PetscErrorCode ISStrideSetStride(IS is,PetscInt n,PetscInt first,PetscInt step)
348: {
351: if (n < 0) SETERRQ1(((PetscObject) is)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Negative length %d not valid", n);
352: PetscUseMethod(is,"ISStrideSetStride_C",(IS,PetscInt,PetscInt,PetscInt),(is,n,first,step));
353: return(0);
354: }
356: EXTERN_C_BEGIN
359: PetscErrorCode ISStrideSetStride_Stride(IS is,PetscInt n,PetscInt first,PetscInt step)
360: {
362: PetscInt min,max;
363: IS_Stride *sub = (IS_Stride*)is->data;
366: sub->n = n;
367: MPI_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,((PetscObject)is)->comm);
368: sub->first = first;
369: sub->step = step;
370: if (step > 0) {min = first; max = first + step*(n-1);}
371: else {max = first; min = first + step*(n-1);}
373: is->min = min;
374: is->max = max;
375: is->data = (void*)sub;
377: if ((!first && step == 1) || (first == max && step == -1 && !min)) {
378: is->isperm = PETSC_TRUE;
379: } else {
380: is->isperm = PETSC_FALSE;
381: }
382: return(0);
383: }
384: EXTERN_C_END
388: /*@
389: ISCreateStride - Creates a data structure for an index set
390: containing a list of evenly spaced integers.
392: Collective on MPI_Comm
394: Input Parameters:
395: + comm - the MPI communicator
396: . n - the length of the locally owned portion of the index set
397: . first - the first element of the locally owned portion of the index set
398: - step - the change to the next index
400: Output Parameter:
401: . is - the new index set
403: Notes:
404: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
405: conceptually the same as MPI_Group operations. The IS are the
406: distributed sets of indices and thus certain operations on them are collective.
408: Level: beginner
410: Concepts: IS^stride
411: Concepts: index sets^stride
412: Concepts: stride^index set
414: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
415: @*/
416: PetscErrorCode ISCreateStride(MPI_Comm comm,PetscInt n,PetscInt first,PetscInt step,IS *is)
417: {
421: ISCreate(comm,is);
422: ISSetType(*is,ISSTRIDE);
423: ISStrideSetStride(*is,n,first,step);
424: return(0);
425: }
427: EXTERN_C_BEGIN
430: PetscErrorCode ISCreate_Stride(IS is)
431: {
433: IS_Stride *sub;
436: PetscMemcpy(is->ops,&myops,sizeof(myops));
437: PetscNewLog(is,IS_Stride,&sub);
438: is->bs = 1;
439: is->data = sub;
440: PetscObjectComposeFunctionDynamic((PetscObject)is,"ISStrideSetStride_C","ISStrideSetStride_Stride",ISStrideSetStride_Stride);
441: return(0);
442: }
443: EXTERN_C_END