Actual source code: stride.c

petsc-master 2019-12-15
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
 61:    the stride width.

 63:    Not Collective

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

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

 72:    Level: intermediate

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


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

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

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

100: PetscErrorCode ISDestroy_Stride(IS is)
101: {

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

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

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

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

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

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

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

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

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

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

178:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
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:   }
218:   return(0);
219: }

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

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

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

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

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

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

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

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

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

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

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

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

282: static PetscErrorCode ISSetBlockSize_Stride(IS is,PetscInt bs)
283: {
284:   IS_Stride     *sub = (IS_Stride*)is->data;

288:   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);
289:   PetscLayoutSetBlockSize(is->map, bs);
290:   return(0);
291: }

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

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


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


335: /*@
336:    ISStrideSetStride - Sets the stride information for a stride index set.

338:    Collective on IS

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

346:    Level: beginner

348: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
349: @*/
350: PetscErrorCode  ISStrideSetStride(IS is,PetscInt n,PetscInt first,PetscInt step)
351: {

355:   if (n < 0) SETERRQ1(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Negative length %D not valid", n);
356:   ISClearInfoCache(is,PETSC_FALSE);
357:   PetscUseMethod(is,"ISStrideSetStride_C",(IS,PetscInt,PetscInt,PetscInt),(is,n,first,step));
358:   return(0);
359: }

361: PetscErrorCode  ISStrideSetStride_Stride(IS is,PetscInt n,PetscInt first,PetscInt step)
362: {
364:   PetscInt       min,max;
365:   IS_Stride      *sub = (IS_Stride*)is->data;
366:   PetscLayout    map;

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

373:   sub->first = first;
374:   sub->step  = step;
375:   if (step > 0) {min = first; max = first + step*(n-1);}
376:   else          {max = first; min = first + step*(n-1);}

378:   is->min  = n > 0 ? min : PETSC_MAX_INT;
379:   is->max  = n > 0 ? max : PETSC_MIN_INT;
380:   is->data = (void*)sub;
381:   return(0);
382: }

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

388:    Collective

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

396:    Output Parameter:
397: .  is - the new index set

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

404:    Level: beginner

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

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

419: PETSC_EXTERN PetscErrorCode ISCreate_Stride(IS is)
420: {
422:   IS_Stride      *sub;

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