Actual source code: stride.c
1: #define PETSCVEC_DLL
2: /*
3: Index sets of evenly space integers, defined by a
4: start, stride and length.
5: */
6: #include src/vec/is/isimpl.h
7: #include petscvec.h
9: typedef struct {
10: PetscInt N,n,first,step;
11: } IS_Stride;
15: PetscErrorCode ISIdentity_Stride(IS is,PetscTruth *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: PetscErrorCode ISDuplicate_Stride(IS is,IS *newIS)
32: {
34: IS_Stride *sub = (IS_Stride*)is->data;
37: ISCreateStride(is->comm,sub->n,sub->first,sub->step,newIS);
38: return(0);
39: }
43: PetscErrorCode ISInvertPermutation_Stride(IS is,PetscInt nlocal,IS *perm)
44: {
45: IS_Stride *isstride = (IS_Stride*)is->data;
49: if (is->isidentity) {
50: ISCreateStride(PETSC_COMM_SELF,isstride->n,0,1,perm);
51: } else {
52: IS tmp;
53: PetscInt *indices,n = isstride->n;
54: ISGetIndices(is,&indices);
55: ISCreateGeneral(is->comm,n,indices,&tmp);
56: ISSetPermutation(tmp);
57: ISRestoreIndices(is,&indices);
58: ISInvertPermutation(tmp,nlocal,perm);
59: ISDestroy(tmp);
60: }
61: return(0);
62: }
63:
66: /*@
67: ISStrideGetInfo - Returns the first index in a stride index set and
68: the stride width.
70: Not Collective
72: Input Parameter:
73: . is - the index set
75: Output Parameters:
76: . first - the first index
77: . step - the stride width
79: Level: intermediate
81: Notes:
82: Returns info on stride index set. This is a pseudo-public function that
83: should not be needed by most users.
85: Concepts: index sets^getting information
86: Concepts: IS^getting information
88: .seealso: ISCreateStride(), ISGetSize()
89: @*/
90: PetscErrorCode ISStrideGetInfo(IS is,PetscInt *first,PetscInt *step)
91: {
92: IS_Stride *sub;
99: sub = (IS_Stride*)is->data;
100: if (is->type != IS_STRIDE) return(0);
101: if (first) *first = sub->first;
102: if (step) *step = sub->step;
103: return(0);
104: }
108: /*@
109: ISStride - Determines if an IS is based on a stride.
111: Not Collective
113: Input Parameter:
114: . is - the index set
116: Output Parameters:
117: . flag - either PETSC_TRUE or PETSC_FALSE
119: Level: intermediate
121: Concepts: index sets^is it stride
122: Concepts: IS^is it stride
124: .seealso: ISCreateStride(), ISGetSize()
125: @*/
126: PetscErrorCode ISStride(IS is,PetscTruth *flag)
127: {
132: if (is->type != IS_STRIDE) *flag = PETSC_FALSE;
133: else *flag = PETSC_TRUE;
135: return(0);
136: }
140: PetscErrorCode ISDestroy_Stride(IS is)
141: {
145: PetscFree(is->data);
146: PetscHeaderDestroy(is);
147: return(0);
148: }
150: /*
151: Returns a legitimate index memory even if
152: the stride index set is empty.
153: */
156: PetscErrorCode ISGetIndices_Stride(IS in,PetscInt **idx)
157: {
158: IS_Stride *sub = (IS_Stride*)in->data;
160: PetscInt i;
163: PetscMalloc(sub->n*sizeof(PetscInt),idx);
164: if (sub->n) {
165: (*idx)[0] = sub->first;
166: for (i=1; i<sub->n; i++) (*idx)[i] = (*idx)[i-1] + sub->step;
167: }
168: return(0);
169: }
173: PetscErrorCode ISRestoreIndices_Stride(IS in,PetscInt **idx)
174: {
178: PetscFree(*idx);
179: return(0);
180: }
184: PetscErrorCode ISGetSize_Stride(IS is,PetscInt *size)
185: {
186: IS_Stride *sub = (IS_Stride *)is->data;
189: *size = sub->N;
190: return(0);
191: }
195: PetscErrorCode ISGetLocalSize_Stride(IS is,PetscInt *size)
196: {
197: IS_Stride *sub = (IS_Stride *)is->data;
200: *size = sub->n;
201: return(0);
202: }
206: PetscErrorCode ISView_Stride(IS is,PetscViewer viewer)
207: {
208: IS_Stride *sub = (IS_Stride *)is->data;
209: PetscInt i,n = sub->n;
210: PetscMPIInt rank,size;
211: PetscTruth iascii;
215: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
216: if (iascii) {
217: MPI_Comm_rank(is->comm,&rank);
218: MPI_Comm_size(is->comm,&size);
219: if (size == 1) {
220: if (is->isperm) {
221: PetscViewerASCIISynchronizedPrintf(viewer,"Index set is permutation\n");
222: }
223: PetscViewerASCIISynchronizedPrintf(viewer,"Number of indices in (stride) set %D\n",n);
224: for (i=0; i<n; i++) {
225: PetscViewerASCIISynchronizedPrintf(viewer,"%D %D\n",i,sub->first + i*sub->step);
226: }
227: } else {
228: if (is->isperm) {
229: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutation\n",rank);
230: }
231: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in (stride) set %D\n",rank,n);
232: for (i=0; i<n; i++) {
233: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,sub->first + i*sub->step);
234: }
235: }
236: PetscViewerFlush(viewer);
237: } else {
238: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for this object",((PetscObject)viewer)->type_name);
239: }
240: return(0);
241: }
242:
245: PetscErrorCode ISSort_Stride(IS is)
246: {
247: IS_Stride *sub = (IS_Stride*)is->data;
250: if (sub->step >= 0) return(0);
251: sub->first += (sub->n - 1)*sub->step;
252: sub->step *= -1;
253: return(0);
254: }
258: PetscErrorCode ISSorted_Stride(IS is,PetscTruth* flg)
259: {
260: IS_Stride *sub = (IS_Stride*)is->data;
263: if (sub->step >= 0) *flg = PETSC_TRUE;
264: else *flg = PETSC_FALSE;
265: return(0);
266: }
268: static struct _ISOps myops = { ISGetSize_Stride,
269: ISGetLocalSize_Stride,
270: ISGetIndices_Stride,
271: ISRestoreIndices_Stride,
272: ISInvertPermutation_Stride,
273: ISSort_Stride,
274: ISSorted_Stride,
275: ISDuplicate_Stride,
276: ISDestroy_Stride,
277: ISView_Stride,
278: ISIdentity_Stride };
282: /*@
283: ISCreateStride - Creates a data structure for an index set
284: containing a list of evenly spaced integers.
286: Collective on MPI_Comm
288: Input Parameters:
289: + comm - the MPI communicator
290: . n - the length of the index set
291: . first - the first element of the index set
292: - step - the change to the next index
294: Output Parameter:
295: . is - the new index set
297: Notes:
298: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
299: conceptually the same as MPI_Group operations. The IS are the
300: distributed sets of indices and thus certain operations on them are collective.
302: Level: beginner
304: Concepts: IS^stride
305: Concepts: index sets^stride
306: Concepts: stride^index set
308: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
309: @*/
310: PetscErrorCode ISCreateStride(MPI_Comm comm,PetscInt n,PetscInt first,PetscInt step,IS *is)
311: {
313: PetscInt min,max;
314: IS Nindex;
315: IS_Stride *sub;
316: PetscTruth flg;
320: *is = PETSC_NULL;
321: if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Number of indices < 0");
322: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
323: VecInitializePackage(PETSC_NULL);
324: #endif
326: PetscHeaderCreate(Nindex,_p_IS,struct _ISOps,IS_COOKIE,IS_STRIDE,"IS",comm,ISDestroy,ISView);
327: PetscLogObjectMemory(Nindex,sizeof(IS_Stride) + sizeof(struct _p_IS));
328: PetscNew(IS_Stride,&sub);
329: sub->n = n;
330: MPI_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,comm);
331: sub->first = first;
332: sub->step = step;
333: if (step > 0) {min = first; max = first + step*(n-1);}
334: else {max = first; min = first + step*(n-1);}
336: Nindex->min = min;
337: Nindex->max = max;
338: Nindex->data = (void*)sub;
339: PetscMemcpy(Nindex->ops,&myops,sizeof(myops));
341: if ((!first && step == 1) || (first == max && step == -1 && !min)) {
342: Nindex->isperm = PETSC_TRUE;
343: } else {
344: Nindex->isperm = PETSC_FALSE;
345: }
346: PetscOptionsHasName(PETSC_NULL,"-is_view",&flg);
347: if (flg) {
348: PetscViewer viewer;
349: PetscViewerASCIIGetStdout(Nindex->comm,&viewer);
350: ISView(Nindex,viewer);
351: }
352: *is = Nindex;
353: return(0);
354: }