Actual source code: stride.c
petsc-3.11.4 2019-09-28
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 <petscvec.h>
8: #include <petscviewer.h>
10: typedef struct {
11: PetscInt N,n,first,step;
12: } IS_Stride;
14: PetscErrorCode ISIdentity_Stride(IS is,PetscBool *ident)
15: {
16: IS_Stride *is_stride = (IS_Stride*)is->data;
19: is->isidentity = PETSC_FALSE;
20: *ident = PETSC_FALSE;
21: if (is_stride->first != 0) return(0);
22: if (is_stride->step != 1) return(0);
23: *ident = PETSC_TRUE;
24: is->isidentity = PETSC_TRUE;
25: return(0);
26: }
28: static PetscErrorCode ISCopy_Stride(IS is,IS isy)
29: {
30: IS_Stride *is_stride = (IS_Stride*)is->data,*isy_stride = (IS_Stride*)isy->data;
34: PetscMemcpy(isy_stride,is_stride,sizeof(IS_Stride));
35: return(0);
36: }
38: PetscErrorCode ISDuplicate_Stride(IS is,IS *newIS)
39: {
41: IS_Stride *sub = (IS_Stride*)is->data;
44: ISCreateStride(PetscObjectComm((PetscObject)is),sub->n,sub->first,sub->step,newIS);
45: return(0);
46: }
48: PetscErrorCode ISInvertPermutation_Stride(IS is,PetscInt nlocal,IS *perm)
49: {
50: IS_Stride *isstride = (IS_Stride*)is->data;
54: if (is->isidentity) {
55: ISCreateStride(PETSC_COMM_SELF,isstride->n,0,1,perm);
56: } else {
57: IS tmp;
58: const PetscInt *indices,n = isstride->n;
59: ISGetIndices(is,&indices);
60: ISCreateGeneral(PetscObjectComm((PetscObject)is),n,indices,PETSC_COPY_VALUES,&tmp);
61: ISSetPermutation(tmp);
62: ISRestoreIndices(is,&indices);
63: ISInvertPermutation(tmp,nlocal,perm);
64: ISDestroy(&tmp);
65: }
66: return(0);
67: }
69: /*@
70: ISStrideGetInfo - Returns the first index in a stride index set and
71: the stride width.
73: Not Collective
75: Input Parameter:
76: . is - the index set
78: Output Parameters:
79: . first - the first index
80: . step - the stride width
82: Level: intermediate
84: Notes:
85: Returns info on stride index set. This is a pseudo-public function that
86: should not be needed by most users.
88: Concepts: index sets^getting information
89: Concepts: IS^getting information
91: .seealso: ISCreateStride(), ISGetSize()
92: @*/
93: PetscErrorCode ISStrideGetInfo(IS is,PetscInt *first,PetscInt *step)
94: {
95: IS_Stride *sub;
96: PetscBool flg;
103: PetscObjectTypeCompare((PetscObject)is,ISSTRIDE,&flg);
104: if (!flg) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_WRONG,"IS must be of type ISSTRIDE");
106: sub = (IS_Stride*)is->data;
107: if (first) *first = sub->first;
108: if (step) *step = sub->step;
109: return(0);
110: }
112: PetscErrorCode ISDestroy_Stride(IS is)
113: {
117: PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",NULL);
118: PetscFree(is->data);
119: return(0);
120: }
122: PetscErrorCode ISToGeneral_Stride(IS inis)
123: {
125: const PetscInt *idx;
126: PetscInt n;
129: ISGetLocalSize(inis,&n);
130: ISGetIndices(inis,&idx);
131: ISSetType(inis,ISGENERAL);
132: ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
133: return(0);
134: }
136: PetscErrorCode ISLocate_Stride(IS is,PetscInt key,PetscInt *location)
137: {
138: IS_Stride *sub = (IS_Stride*)is->data;
139: PetscInt rem, step;
142: *location = -1;
143: step = sub->step;
144: key -= sub->first;
145: rem = key / step;
146: if ((rem < sub->n) && !(key % step)) {
147: *location = rem;
148: }
149: return(0);
150: }
152: /*
153: Returns a legitimate index memory even if
154: the stride index set is empty.
155: */
156: PetscErrorCode ISGetIndices_Stride(IS in,const PetscInt *idx[])
157: {
158: IS_Stride *sub = (IS_Stride*)in->data;
160: PetscInt i,**dx = (PetscInt**)idx;
163: PetscMalloc1(sub->n,(PetscInt**)idx);
164: if (sub->n) {
165: (*dx)[0] = sub->first;
166: for (i=1; i<sub->n; i++) (*dx)[i] = (*dx)[i-1] + sub->step;
167: }
168: return(0);
169: }
171: PetscErrorCode ISRestoreIndices_Stride(IS in,const PetscInt *idx[])
172: {
176: PetscFree(*(void**)idx);
177: return(0);
178: }
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: }
189: PetscErrorCode ISGetLocalSize_Stride(IS is,PetscInt *size)
190: {
191: IS_Stride *sub = (IS_Stride*)is->data;
194: *size = sub->n;
195: return(0);
196: }
198: PetscErrorCode ISView_Stride(IS is,PetscViewer viewer)
199: {
200: IS_Stride *sub = (IS_Stride*)is->data;
201: PetscInt i,n = sub->n;
202: PetscMPIInt rank,size;
203: PetscBool iascii;
204: PetscViewerFormat fmt;
205: PetscErrorCode ierr;
208: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
209: if (iascii) {
210: PetscBool matl;
212: MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
213: MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
214: PetscViewerGetFormat(viewer,&fmt);
215: matl = (PetscBool)(fmt == PETSC_VIEWER_ASCII_MATLAB);
216: if (size == 1) {
217: if (matl) {
218: const char* name;
220: PetscObjectGetName((PetscObject)is,&name);
221: PetscViewerASCIIPrintf(viewer,"%s = [%D : %D : %D];\n",name,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
222: } else {
223: if (is->isperm) {
224: PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");
225: }
226: PetscViewerASCIIPrintf(viewer,"Number of indices in (stride) set %D\n",n);
227: for (i=0; i<n; i++) {
228: PetscViewerASCIIPrintf(viewer,"%D %D\n",i,sub->first + i*sub->step);
229: }
230: }
231: PetscViewerFlush(viewer);
232: } else {
233: PetscViewerASCIIPushSynchronized(viewer);
234: if (matl) {
235: const char* name;
237: PetscObjectGetName((PetscObject)is,&name);
238: PetscViewerASCIISynchronizedPrintf(viewer,"%s_%d = [%D : %D : %D];\n",name,rank,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
239: } else {
240: if (is->isperm) {
241: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutation\n",rank);
242: }
243: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in (stride) set %D\n",rank,n);
244: for (i=0; i<n; i++) {
245: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,sub->first + i*sub->step);
246: }
247: }
248: PetscViewerFlush(viewer);
249: PetscViewerASCIIPopSynchronized(viewer);
250: }
251: }
252: return(0);
253: }
255: PetscErrorCode ISSort_Stride(IS is)
256: {
257: IS_Stride *sub = (IS_Stride*)is->data;
260: if (sub->step >= 0) return(0);
261: sub->first += (sub->n - 1)*sub->step;
262: sub->step *= -1;
263: return(0);
264: }
266: PetscErrorCode ISSorted_Stride(IS is,PetscBool * flg)
267: {
268: IS_Stride *sub = (IS_Stride*)is->data;
271: if (sub->step >= 0) *flg = PETSC_TRUE;
272: else *flg = PETSC_FALSE;
273: return(0);
274: }
276: static PetscErrorCode ISOnComm_Stride(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
277: {
279: IS_Stride *sub = (IS_Stride*)is->data;
282: ISCreateStride(comm,sub->n,sub->first,sub->step,newis);
283: return(0);
284: }
286: static PetscErrorCode ISSetBlockSize_Stride(IS is,PetscInt bs)
287: {
288: IS_Stride *sub = (IS_Stride*)is->data;
292: 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);
293: PetscLayoutSetBlockSize(is->map, bs);
294: return(0);
295: }
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,
331: ISLocate_Stride};
334: /*@
335: ISStrideSetStride - Sets the stride information for a stride index set.
337: Collective on IS
339: Input Parameters:
340: + is - the index set
341: . n - the length of the locally owned portion of the index set
342: . first - the first element of the locally owned portion of the index set
343: - step - the change to the next index
345: Level: beginner
347: Concepts: IS^stride
348: Concepts: index sets^stride
349: Concepts: stride^index set
351: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
352: @*/
353: PetscErrorCode ISStrideSetStride(IS is,PetscInt n,PetscInt first,PetscInt step)
354: {
358: if (n < 0) SETERRQ1(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Negative length %d not valid", n);
359: PetscUseMethod(is,"ISStrideSetStride_C",(IS,PetscInt,PetscInt,PetscInt),(is,n,first,step));
360: return(0);
361: }
363: PetscErrorCode ISStrideSetStride_Stride(IS is,PetscInt n,PetscInt first,PetscInt step)
364: {
366: PetscInt min,max;
367: IS_Stride *sub = (IS_Stride*)is->data;
370: sub->n = n;
371: MPIU_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)is));
372: sub->first = first;
373: sub->step = step;
374: if (step > 0) {min = first; max = first + step*(n-1);}
375: else {max = first; min = first + step*(n-1);}
377: is->min = n > 0 ? min : PETSC_MAX_INT;
378: is->max = n > 0 ? max : PETSC_MIN_INT;
379: is->data = (void*)sub;
381: if ((!first && step == 1) || (first == max && step == -1 && !min)) is->isperm = PETSC_TRUE;
382: else is->isperm = PETSC_FALSE;
383: is->isidentity = PETSC_FALSE;
384: return(0);
385: }
387: /*@
388: ISCreateStride - Creates a data structure for an index set
389: containing a list of evenly spaced integers.
391: Collective on MPI_Comm
393: Input Parameters:
394: + comm - the MPI communicator
395: . n - the length of the locally owned portion of the index set
396: . first - the first element of the locally owned portion of the index set
397: - step - the change to the next index
399: Output Parameter:
400: . is - the new index set
402: Notes:
403: When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
404: conceptually the same as MPI_Group operations. The IS are the
405: distributed sets of indices and thus certain operations on them are collective.
407: Level: beginner
409: Concepts: IS^stride
410: Concepts: index sets^stride
411: Concepts: stride^index set
413: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
414: @*/
415: PetscErrorCode ISCreateStride(MPI_Comm comm,PetscInt n,PetscInt first,PetscInt step,IS *is)
416: {
420: ISCreate(comm,is);
421: ISSetType(*is,ISSTRIDE);
422: ISStrideSetStride(*is,n,first,step);
423: return(0);
424: }
426: PETSC_EXTERN PetscErrorCode ISCreate_Stride(IS is)
427: {
429: IS_Stride *sub;
432: PetscNewLog(is,&sub);
433: is->data = (void *) sub;
434: PetscMemcpy(is->ops,&myops,sizeof(myops));
435: PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",ISStrideSetStride_Stride);
436: return(0);
437: }