Actual source code: stride.c
1: /*$Id: stride.c,v 1.104 2001/06/21 21:15:53 bsmith Exp $*/
2: /*
3: Index sets of evenly space integers, defined by a
4: start, stride and length.
5: */
6: #include src/vec/is/isimpl.h
8: EXTERN int VecInitializePackage(char *);
10: typedef struct {
11: int N,n,first,step;
12: } IS_Stride;
14: #undef __FUNCT__
16: int ISIdentity_Stride(IS is,PetscTruth *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: }
30: #undef __FUNCT__
32: int ISDuplicate_Stride(IS is,IS *newIS)
33: {
34: int ierr;
35: IS_Stride *sub = (IS_Stride*)is->data;
38: ISCreateStride(is->comm,sub->n,sub->first,sub->step,newIS);
39: return(0);
40: }
42: #undef __FUNCT__
44: int ISInvertPermutation_Stride(IS is,int nlocal,IS *perm)
45: {
46: IS_Stride *isstride = (IS_Stride*)is->data;
47: int ierr;
50: if (is->isidentity) {
51: ISCreateStride(PETSC_COMM_SELF,isstride->n,0,1,perm);
52: } else {
53: IS tmp;
54: int *indices,n = isstride->n;
55: ISGetIndices(is,&indices);
56: ISCreateGeneral(is->comm,n,indices,&tmp);
57: ISRestoreIndices(is,&indices);
58: ISInvertPermutation(tmp,nlocal,perm);
59: ISDestroy(tmp);
60: }
61: return(0);
62: }
63:
64: #undef __FUNCT__
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: int ISStrideGetInfo(IS is,int *first,int *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: }
106: #undef __FUNCT__
108: /*@C
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: int ISStride(IS is,PetscTruth *flag)
127: {
132: if (is->type != IS_STRIDE) *flag = PETSC_FALSE;
133: else *flag = PETSC_TRUE;
135: return(0);
136: }
138: #undef __FUNCT__
140: int ISDestroy_Stride(IS is)
141: {
145: PetscFree(is->data);
146: PetscLogObjectDestroy(is);
147: PetscHeaderDestroy(is); return(0);
148: }
150: /*
151: Returns a legitimate index memory even if
152: the stride index set is empty.
153: */
154: #undef __FUNCT__
156: int ISGetIndices_Stride(IS in,int **idx)
157: {
158: IS_Stride *sub = (IS_Stride*)in->data;
159: int i,ierr;
162: ierr = PetscMalloc((sub->n+1)*sizeof(int),idx);
163: (*idx)[0] = sub->first;
164: for (i=1; i<sub->n; i++) (*idx)[i] = (*idx)[i-1] + sub->step;
165: return(0);
166: }
168: #undef __FUNCT__
170: int ISRestoreIndices_Stride(IS in,int **idx)
171: {
175: PetscFree(*idx);
176: return(0);
177: }
179: #undef __FUNCT__
181: int ISGetSize_Stride(IS is,int *size)
182: {
183: IS_Stride *sub = (IS_Stride *)is->data;
186: *size = sub->N;
187: return(0);
188: }
190: #undef __FUNCT__
192: int ISGetLocalSize_Stride(IS is,int *size)
193: {
194: IS_Stride *sub = (IS_Stride *)is->data;
197: *size = sub->n;
198: return(0);
199: }
201: #undef __FUNCT__
203: int ISView_Stride(IS is,PetscViewer viewer)
204: {
205: IS_Stride *sub = (IS_Stride *)is->data;
206: int i,n = sub->n,ierr,rank,size;
207: PetscTruth isascii;
210: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
211: if (isascii) {
212: MPI_Comm_rank(is->comm,&rank);
213: MPI_Comm_size(is->comm,&size);
214: if (size == 1) {
215: if (is->isperm) {
216: PetscViewerASCIISynchronizedPrintf(viewer,"Index set is permutationn");
217: }
218: PetscViewerASCIISynchronizedPrintf(viewer,"Number of indices in (stride) set %dn",n);
219: for (i=0; i<n; i++) {
220: PetscViewerASCIISynchronizedPrintf(viewer,"%d %dn",i,sub->first + i*sub->step);
221: }
222: } else {
223: if (is->isperm) {
224: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutationn",rank);
225: }
226: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in (stride) set %dn",rank,n);
227: for (i=0; i<n; i++) {
228: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %d %dn",rank,i,sub->first + i*sub->step);
229: }
230: }
231: PetscViewerFlush(viewer);
232: } else {
233: SETERRQ1(1,"Viewer type %s not supported for this object",((PetscObject)viewer)->type_name);
234: }
235: return(0);
236: }
237:
238: #undef __FUNCT__
240: int ISSort_Stride(IS is)
241: {
242: IS_Stride *sub = (IS_Stride*)is->data;
245: if (sub->step >= 0) return(0);
246: sub->first += (sub->n - 1)*sub->step;
247: sub->step *= -1;
248: return(0);
249: }
251: #undef __FUNCT__
253: int ISSorted_Stride(IS is,PetscTruth* flg)
254: {
255: IS_Stride *sub = (IS_Stride*)is->data;
258: if (sub->step >= 0) *flg = PETSC_TRUE;
259: else *flg = PETSC_FALSE;
260: return(0);
261: }
263: static struct _ISOps myops = { ISGetSize_Stride,
264: ISGetLocalSize_Stride,
265: ISGetIndices_Stride,
266: ISRestoreIndices_Stride,
267: ISInvertPermutation_Stride,
268: ISSort_Stride,
269: ISSorted_Stride,
270: ISDuplicate_Stride,
271: ISDestroy_Stride,
272: ISView_Stride,
273: ISIdentity_Stride };
275: #undef __FUNCT__
277: /*@C
278: ISCreateStride - Creates a data structure for an index set
279: containing a list of evenly spaced integers.
281: Collective on MPI_Comm
283: Input Parameters:
284: + comm - the MPI communicator
285: . n - the length of the index set
286: . first - the first element of the index set
287: - step - the change to the next index
289: Output Parameter:
290: . is - the new index set
292: Notes:
293: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
294: conceptually the same as MPI_Group operations. The IS are the
295: distributed sets of indices and thus certain operations on them are collective.
297: Level: beginner
299: Concepts: IS^stride
300: Concepts: index sets^stride
301: Concepts: stride^index set
303: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
304: @*/
305: int ISCreateStride(MPI_Comm comm,int n,int first,int step,IS *is)
306: {
307: int min,max,ierr;
308: IS Nindex;
309: IS_Stride *sub;
310: PetscTruth flg;
314: *is = PETSC_NULL;
315: if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Number of indices < 0");
316: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
317: VecInitializePackage(PETSC_NULL);
318: #endif
320: PetscHeaderCreate(Nindex,_p_IS,struct _ISOps,IS_COOKIE,IS_STRIDE,"IS",comm,ISDestroy,ISView);
321: PetscLogObjectCreate(Nindex);
322: PetscLogObjectMemory(Nindex,sizeof(IS_Stride) + sizeof(struct _p_IS));
323: ierr = PetscNew(IS_Stride,&sub);
324: sub->n = n;
325: MPI_Allreduce(&n,&sub->N,1,MPI_INT,MPI_SUM,comm);
326: sub->first = first;
327: sub->step = step;
328: if (step > 0) {min = first; max = first + step*(n-1);}
329: else {max = first; min = first + step*(n-1);}
331: Nindex->min = min;
332: Nindex->max = max;
333: Nindex->data = (void*)sub;
334: PetscMemcpy(Nindex->ops,&myops,sizeof(myops));
336: if ((first == 0 && step == 1) || (first == max && step == -1 && min == 0)) {
337: Nindex->isperm = PETSC_TRUE;
338: } else {
339: Nindex->isperm = PETSC_FALSE;
340: }
341: PetscOptionsHasName(PETSC_NULL,"-is_view",&flg);
342: if (flg) {
343: ISView(Nindex,PETSC_VIEWER_STDOUT_(Nindex->comm));
344: }
345: *is = Nindex;
346: return(0);
347: }