Actual source code: stride.c

petsc-master 2020-12-02
Report Typos and Errors

  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 <petscviewer.h>

  9: typedef struct {
 10:   PetscInt first,step;
 11: } IS_Stride;

 13: static PetscErrorCode ISCopy_Stride(IS is,IS isy)
 14: {
 15:   IS_Stride      *is_stride = (IS_Stride*)is->data,*isy_stride = (IS_Stride*)isy->data;

 19:   PetscMemcpy(isy_stride,is_stride,sizeof(IS_Stride));
 20:   return(0);
 21: }

 23: PetscErrorCode ISDuplicate_Stride(IS is,IS *newIS)
 24: {
 26:   IS_Stride      *sub = (IS_Stride*)is->data;

 29:   ISCreateStride(PetscObjectComm((PetscObject)is),is->map->n,sub->first,sub->step,newIS);
 30:   return(0);
 31: }

 33: PetscErrorCode ISInvertPermutation_Stride(IS is,PetscInt nlocal,IS *perm)
 34: {
 35:   PetscBool      isident;

 39:   ISGetInfo(is,IS_IDENTITY,IS_GLOBAL,PETSC_TRUE,&isident);
 40:   if (isident) {
 41:     PetscInt rStart, rEnd;

 43:     PetscLayoutGetRange(is->map, &rStart, &rEnd);
 44:     ISCreateStride(PETSC_COMM_SELF,PetscMax(rEnd - rStart, 0),rStart,1,perm);
 45:   } else {
 46:     IS             tmp;
 47:     const PetscInt *indices,n = is->map->n;

 49:     ISGetIndices(is,&indices);
 50:     ISCreateGeneral(PetscObjectComm((PetscObject)is),n,indices,PETSC_COPY_VALUES,&tmp);
 51:     ISSetPermutation(tmp);
 52:     ISRestoreIndices(is,&indices);
 53:     ISInvertPermutation(tmp,nlocal,perm);
 54:     ISDestroy(&tmp);
 55:   }
 56:   return(0);
 57: }

 59: /*@
 60:    ISStrideGetInfo - Returns the first index in a stride index set and the stride width.

 62:    Not Collective

 64:    Input Parameter:
 65: .  is - the index set

 67:    Output Parameters:
 68: +  first - the first index
 69: -  step - the stride width

 71:    Level: intermediate

 73:    Notes:
 74:    Returns info on stride index set. This is a pseudo-public function that
 75:    should not be needed by most users.


 78: .seealso: ISCreateStride(), ISGetSize(), ISSTRIDE
 79: @*/
 80: PetscErrorCode  ISStrideGetInfo(IS is,PetscInt *first,PetscInt *step)
 81: {
 82:   IS_Stride      *sub;
 83:   PetscBool      flg;

 90:   PetscObjectTypeCompare((PetscObject)is,ISSTRIDE,&flg);
 91:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_WRONG,"IS must be of type ISSTRIDE");

 93:   sub = (IS_Stride*)is->data;
 94:   if (first) *first = sub->first;
 95:   if (step)  *step  = sub->step;
 96:   return(0);
 97: }

 99: PetscErrorCode ISDestroy_Stride(IS is)
100: {

104:   PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",NULL);
105:   PetscFree(is->data);
106:   return(0);
107: }

109: PetscErrorCode  ISToGeneral_Stride(IS inis)
110: {
112:   const PetscInt *idx;
113:   PetscInt       n;

116:   ISGetLocalSize(inis,&n);
117:   ISGetIndices(inis,&idx);
118:   ISSetType(inis,ISGENERAL);
119:   ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
120:   return(0);
121: }

123: PetscErrorCode ISLocate_Stride(IS is,PetscInt key,PetscInt *location)
124: {
125:   IS_Stride      *sub = (IS_Stride*)is->data;
126:   PetscInt       rem, step;

129:   *location = -1;
130:   step      = sub->step;
131:   key      -= sub->first;
132:   rem       = key / step;
133:   if ((rem < is->map->n) && !(key % step)) {
134:     *location = rem;
135:   }
136:   return(0);
137: }

139: /*
140:      Returns a legitimate index memory even if
141:    the stride index set is empty.
142: */
143: PetscErrorCode ISGetIndices_Stride(IS is,const PetscInt *idx[])
144: {
145:   IS_Stride      *sub = (IS_Stride*)is->data;
147:   PetscInt       i,**dx = (PetscInt**)idx;

150:   PetscMalloc1(is->map->n,(PetscInt**)idx);
151:   if (is->map->n) {
152:     (*dx)[0] = sub->first;
153:     for (i=1; i<is->map->n; i++) (*dx)[i] = (*dx)[i-1] + sub->step;
154:   }
155:   return(0);
156: }

158: PetscErrorCode ISRestoreIndices_Stride(IS in,const PetscInt *idx[])
159: {

163:   PetscFree(*(void**)idx);
164:   return(0);
165: }

167: PetscErrorCode ISView_Stride(IS is,PetscViewer viewer)
168: {
169:   IS_Stride         *sub = (IS_Stride*)is->data;
170:   PetscInt          i,n = is->map->n;
171:   PetscMPIInt       rank,size;
172:   PetscBool         iascii,ibinary;
173:   PetscViewerFormat fmt;
174:   PetscErrorCode    ierr;

177:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
178:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&ibinary);
179:   if (iascii) {
180:     PetscBool matl, isperm;

182:     MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
183:     MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
184:     PetscViewerGetFormat(viewer,&fmt);
185:     matl = (PetscBool)(fmt == PETSC_VIEWER_ASCII_MATLAB);
186:     ISGetInfo(is,IS_PERMUTATION,IS_GLOBAL,PETSC_FALSE,&isperm);
187:     if (isperm && !matl) {PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");}
188:     if (size == 1) {
189:       if (matl) {
190:         const char* name;

192:         PetscObjectGetName((PetscObject)is,&name);
193:         PetscViewerASCIIPrintf(viewer,"%s = [%D : %D : %D];\n",name,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
194:       } else {
195:         PetscViewerASCIIPrintf(viewer,"Number of indices in (stride) set %D\n",n);
196:         for (i=0; i<n; i++) {
197:           PetscViewerASCIIPrintf(viewer,"%D %D\n",i,sub->first + i*sub->step);
198:         }
199:       }
200:       PetscViewerFlush(viewer);
201:     } else {
202:       PetscViewerASCIIPushSynchronized(viewer);
203:       if (matl) {
204:         const char* name;

206:         PetscObjectGetName((PetscObject)is,&name);
207:         PetscViewerASCIISynchronizedPrintf(viewer,"%s_%d = [%D : %D : %D];\n",name,rank,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
208:       } else {
209:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in (stride) set %D\n",rank,n);
210:         for (i=0; i<n; i++) {
211:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,sub->first + i*sub->step);
212:         }
213:       }
214:       PetscViewerFlush(viewer);
215:       PetscViewerASCIIPopSynchronized(viewer);
216:     }
217:   } else if (ibinary) {
218:     ISView_Binary(is,viewer);
219:   }
220:   return(0);
221: }

223: PetscErrorCode ISSort_Stride(IS is)
224: {
225:   IS_Stride *sub = (IS_Stride*)is->data;

228:   if (sub->step >= 0) return(0);
229:   sub->first += (is->map->n - 1)*sub->step;
230:   sub->step  *= -1;
231:   return(0);
232: }

234: PetscErrorCode ISSorted_Stride(IS is,PetscBool * flg)
235: {
236:   IS_Stride *sub = (IS_Stride*)is->data;

239:   if (sub->step >= 0) *flg = PETSC_TRUE;
240:   else *flg = PETSC_FALSE;
241:   return(0);
242: }

244: static PetscErrorCode ISUniqueLocal_Stride(IS is, PetscBool *flg)
245: {
246:   IS_Stride *sub = (IS_Stride*)is->data;

249:   if (!(is->map->n) || sub->step != 0) *flg = PETSC_TRUE;
250:   else *flg = PETSC_FALSE;
251:   return(0);
252: }

254: static PetscErrorCode ISPermutationLocal_Stride(IS is, PetscBool *flg)
255: {
256:   IS_Stride *sub = (IS_Stride*)is->data;

259:   if (!(is->map->n) || (PetscAbsInt(sub->step) == 1 && is->min == 0)) *flg = PETSC_TRUE;
260:   else *flg = PETSC_FALSE;
261:   return(0);
262: }

264: static PetscErrorCode ISIntervalLocal_Stride(IS is, PetscBool *flg)
265: {
266:   IS_Stride *sub = (IS_Stride*)is->data;

269:   if (!(is->map->n) || sub->step == 1) *flg = PETSC_TRUE;
270:   else *flg = PETSC_FALSE;
271:   return(0);
272: }

274: static PetscErrorCode ISOnComm_Stride(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
275: {
277:   IS_Stride      *sub = (IS_Stride*)is->data;

280:   ISCreateStride(comm,is->map->n,sub->first,sub->step,newis);
281:   return(0);
282: }

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: }

295: static PetscErrorCode ISContiguousLocal_Stride(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
296: {
297:   IS_Stride *sub = (IS_Stride*)is->data;

300:   if (sub->step == 1 && sub->first >= gstart && sub->first+is->map->n <= gend) {
301:     *start  = sub->first - gstart;
302:     *contig = PETSC_TRUE;
303:   } else {
304:     *start  = -1;
305:     *contig = PETSC_FALSE;
306:   }
307:   return(0);
308: }


311: static struct _ISOps myops = { ISGetIndices_Stride,
312:                                ISRestoreIndices_Stride,
313:                                ISInvertPermutation_Stride,
314:                                ISSort_Stride,
315:                                ISSort_Stride,
316:                                ISSorted_Stride,
317:                                ISDuplicate_Stride,
318:                                ISDestroy_Stride,
319:                                ISView_Stride,
320:                                ISLoad_Default,
321:                                ISCopy_Stride,
322:                                ISToGeneral_Stride,
323:                                ISOnComm_Stride,
324:                                ISSetBlockSize_Stride,
325:                                ISContiguousLocal_Stride,
326:                                ISLocate_Stride,
327:                                ISSorted_Stride,
328:                                NULL,
329:                                ISUniqueLocal_Stride,
330:                                NULL,
331:                                ISPermutationLocal_Stride,
332:                                NULL,
333:                                ISIntervalLocal_Stride,
334:                                NULL};


337: /*@
338:    ISStrideSetStride - Sets the stride information for a stride index set.

340:    Collective on IS

342:    Input Parameters:
343: +  is - the index set
344: .  n - the length of the locally owned portion of the index set
345: .  first - the first element of the locally owned portion of the index set
346: -  step - the change to the next index

348:    Level: beginner

350: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather(), ISSTRIDE, ISCreateStride(), ISStrideGetInfo()
351: @*/
352: PetscErrorCode  ISStrideSetStride(IS is,PetscInt n,PetscInt first,PetscInt step)
353: {

357:   if (n < 0) SETERRQ1(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Negative length %D not valid", n);
358:   ISClearInfoCache(is,PETSC_FALSE);
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;
368:   PetscLayout    map;

371:   PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is),n,is->map->N,is->map->bs,&map);
372:   PetscLayoutDestroy(&is->map);
373:   is->map = map;

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 : PETSC_MAX_INT;
381:   is->max  = n > 0 ? max : PETSC_MIN_INT;
382:   is->data = (void*)sub;
383:   return(0);
384: }

386: /*@
387:    ISCreateStride - Creates a data structure for an index set
388:    containing a list of evenly spaced integers.

390:    Collective

392:    Input Parameters:
393: +  comm - the MPI communicator
394: .  n - the length of the locally owned portion of the index set
395: .  first - the first element of the locally owned portion of the index set
396: -  step - the change to the next index

398:    Output Parameter:
399: .  is - the new index set

401:    Notes:
402:    When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
403:    conceptually the same as MPI_Group operations. The IS are the
404:    distributed sets of indices and thus certain operations on them are collective.

406:    Level: beginner

408: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather(), ISSTRIDE
409: @*/
410: PetscErrorCode  ISCreateStride(MPI_Comm comm,PetscInt n,PetscInt first,PetscInt step,IS *is)
411: {

415:   ISCreate(comm,is);
416:   ISSetType(*is,ISSTRIDE);
417:   ISStrideSetStride(*is,n,first,step);
418:   return(0);
419: }

421: PETSC_EXTERN PetscErrorCode ISCreate_Stride(IS is)
422: {
424:   IS_Stride      *sub;

427:   PetscNewLog(is,&sub);
428:   is->data = (void *) sub;
429:   PetscMemcpy(is->ops,&myops,sizeof(myops));
430:   PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",ISStrideSetStride_Stride);
431:   return(0);
432: }