Actual source code: stride.c

petsc-master 2019-05-22
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 N,n,first,step;
 11: } IS_Stride;

 13: PetscErrorCode ISIdentity_Stride(IS is,PetscBool  *ident)
 14: {
 15:   IS_Stride *is_stride = (IS_Stride*)is->data;

 18:   is->isidentity = PETSC_FALSE;
 19:   *ident         = PETSC_FALSE;
 20:   if (is_stride->first != 0) return(0);
 21:   if (is_stride->step  != 1) return(0);
 22:   *ident         = PETSC_TRUE;
 23:   is->isidentity = PETSC_TRUE;
 24:   return(0);
 25: }

 27: static PetscErrorCode ISCopy_Stride(IS is,IS isy)
 28: {
 29:   IS_Stride      *is_stride = (IS_Stride*)is->data,*isy_stride = (IS_Stride*)isy->data;

 33:   PetscMemcpy(isy_stride,is_stride,sizeof(IS_Stride));
 34:   return(0);
 35: }

 37: PetscErrorCode ISDuplicate_Stride(IS is,IS *newIS)
 38: {
 40:   IS_Stride      *sub = (IS_Stride*)is->data;

 43:   ISCreateStride(PetscObjectComm((PetscObject)is),sub->n,sub->first,sub->step,newIS);
 44:   return(0);
 45: }

 47: PetscErrorCode ISInvertPermutation_Stride(IS is,PetscInt nlocal,IS *perm)
 48: {
 49:   IS_Stride      *isstride = (IS_Stride*)is->data;

 53:   if (is->isidentity) {
 54:     ISCreateStride(PETSC_COMM_SELF,isstride->n,0,1,perm);
 55:   } else {
 56:     IS             tmp;
 57:     const PetscInt *indices,n = isstride->n;
 58:     ISGetIndices(is,&indices);
 59:     ISCreateGeneral(PetscObjectComm((PetscObject)is),n,indices,PETSC_COPY_VALUES,&tmp);
 60:     ISSetPermutation(tmp);
 61:     ISRestoreIndices(is,&indices);
 62:     ISInvertPermutation(tmp,nlocal,perm);
 63:     ISDestroy(&tmp);
 64:   }
 65:   return(0);
 66: }

 68: /*@
 69:    ISStrideGetInfo - Returns the first index in a stride index set and
 70:    the stride width.

 72:    Not Collective

 74:    Input Parameter:
 75: .  is - the index set

 77:    Output Parameters:
 78: .  first - the first index
 79: .  step - the stride width

 81:    Level: intermediate

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

 87:    Concepts: index sets^getting information
 88:    Concepts: IS^getting information

 90: .seealso: ISCreateStride(), ISGetSize()
 91: @*/
 92: PetscErrorCode  ISStrideGetInfo(IS is,PetscInt *first,PetscInt *step)
 93: {
 94:   IS_Stride      *sub;
 95:   PetscBool      flg;

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

105:   sub = (IS_Stride*)is->data;
106:   if (first) *first = sub->first;
107:   if (step)  *step  = sub->step;
108:   return(0);
109: }

111: PetscErrorCode ISDestroy_Stride(IS is)
112: {

116:   PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",NULL);
117:   PetscFree(is->data);
118:   return(0);
119: }

121: PetscErrorCode  ISToGeneral_Stride(IS inis)
122: {
124:   const PetscInt *idx;
125:   PetscInt       n;

128:   ISGetLocalSize(inis,&n);
129:   ISGetIndices(inis,&idx);
130:   ISSetType(inis,ISGENERAL);
131:   ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
132:   return(0);
133: }

135: PetscErrorCode ISLocate_Stride(IS is,PetscInt key,PetscInt *location)
136: {
137:   IS_Stride      *sub = (IS_Stride*)is->data;
138:   PetscInt       rem, step;

141:   *location = -1;
142:   step      = sub->step;
143:   key      -= sub->first;
144:   rem       = key / step;
145:   if ((rem < sub->n) && !(key % step)) {
146:     *location = rem;
147:   }
148:   return(0);
149: }

151: /*
152:      Returns a legitimate index memory even if
153:    the stride index set is empty.
154: */
155: PetscErrorCode ISGetIndices_Stride(IS in,const PetscInt *idx[])
156: {
157:   IS_Stride      *sub = (IS_Stride*)in->data;
159:   PetscInt       i,**dx = (PetscInt**)idx;

162:   PetscMalloc1(sub->n,(PetscInt**)idx);
163:   if (sub->n) {
164:     (*dx)[0] = sub->first;
165:     for (i=1; i<sub->n; i++) (*dx)[i] = (*dx)[i-1] + sub->step;
166:   }
167:   return(0);
168: }

170: PetscErrorCode ISRestoreIndices_Stride(IS in,const PetscInt *idx[])
171: {

175:   PetscFree(*(void**)idx);
176:   return(0);
177: }

179: PetscErrorCode ISGetSize_Stride(IS is,PetscInt *size)
180: {
181:   IS_Stride *sub = (IS_Stride*)is->data;

184:   *size = sub->N;
185:   return(0);
186: }

188: PetscErrorCode ISGetLocalSize_Stride(IS is,PetscInt *size)
189: {
190:   IS_Stride *sub = (IS_Stride*)is->data;

193:   *size = sub->n;
194:   return(0);
195: }

197: PetscErrorCode ISView_Stride(IS is,PetscViewer viewer)
198: {
199:   IS_Stride         *sub = (IS_Stride*)is->data;
200:   PetscInt          i,n = sub->n;
201:   PetscMPIInt       rank,size;
202:   PetscBool         iascii;
203:   PetscViewerFormat fmt;
204:   PetscErrorCode    ierr;

207:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
208:   if (iascii) {
209:     PetscBool matl;

211:     MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
212:     MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
213:     PetscViewerGetFormat(viewer,&fmt);
214:     matl = (PetscBool)(fmt == PETSC_VIEWER_ASCII_MATLAB);
215:     if (size == 1) {
216:       if (matl) {
217:         const char* name;

219:         PetscObjectGetName((PetscObject)is,&name);
220:         PetscViewerASCIIPrintf(viewer,"%s = [%D : %D : %D];\n",name,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
221:       } else {
222:         if (is->isperm) {
223:           PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");
224:         }
225:         PetscViewerASCIIPrintf(viewer,"Number of indices in (stride) set %D\n",n);
226:         for (i=0; i<n; i++) {
227:           PetscViewerASCIIPrintf(viewer,"%D %D\n",i,sub->first + i*sub->step);
228:         }
229:       }
230:       PetscViewerFlush(viewer);
231:     } else {
232:       PetscViewerASCIIPushSynchronized(viewer);
233:       if (matl) {
234:         const char* name;

236:         PetscObjectGetName((PetscObject)is,&name);
237:         PetscViewerASCIISynchronizedPrintf(viewer,"%s_%d = [%D : %D : %D];\n",name,rank,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
238:       } else {
239:         if (is->isperm) {
240:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutation\n",rank);
241:         }
242:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in (stride) set %D\n",rank,n);
243:         for (i=0; i<n; i++) {
244:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,sub->first + i*sub->step);
245:         }
246:       }
247:       PetscViewerFlush(viewer);
248:       PetscViewerASCIIPopSynchronized(viewer);
249:     }
250:   }
251:   return(0);
252: }

254: PetscErrorCode ISSort_Stride(IS is)
255: {
256:   IS_Stride *sub = (IS_Stride*)is->data;

259:   if (sub->step >= 0) return(0);
260:   sub->first += (sub->n - 1)*sub->step;
261:   sub->step  *= -1;
262:   return(0);
263: }

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

270:   if (sub->step >= 0) *flg = PETSC_TRUE;
271:   else *flg = PETSC_FALSE;
272:   return(0);
273: }

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

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

285: static PetscErrorCode ISSetBlockSize_Stride(IS is,PetscInt bs)
286: {
287:   IS_Stride     *sub = (IS_Stride*)is->data;

291:   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);
292:   PetscLayoutSetBlockSize(is->map, bs);
293:   return(0);
294: }

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

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


312: static struct _ISOps myops = { ISGetSize_Stride,
313:                                ISGetLocalSize_Stride,
314:                                ISGetIndices_Stride,
315:                                ISRestoreIndices_Stride,
316:                                ISInvertPermutation_Stride,
317:                                ISSort_Stride,
318:                                ISSort_Stride,
319:                                ISSorted_Stride,
320:                                ISDuplicate_Stride,
321:                                ISDestroy_Stride,
322:                                ISView_Stride,
323:                                ISLoad_Default,
324:                                ISIdentity_Stride,
325:                                ISCopy_Stride,
326:                                ISToGeneral_Stride,
327:                                ISOnComm_Stride,
328:                                ISSetBlockSize_Stride,
329:                                ISContiguousLocal_Stride,
330:                                ISLocate_Stride};


333: /*@
334:    ISStrideSetStride - Sets the stride information for a stride index set.

336:    Collective on IS

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

344:    Level: beginner

346:   Concepts: IS^stride
347:   Concepts: index sets^stride
348:   Concepts: stride^index set

350: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
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:   PetscUseMethod(is,"ISStrideSetStride_C",(IS,PetscInt,PetscInt,PetscInt),(is,n,first,step));
359:   return(0);
360: }

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

369:   sub->n     = n;
370:   MPIU_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)is));
371:   sub->first = first;
372:   sub->step  = step;
373:   if (step > 0) {min = first; max = first + step*(n-1);}
374:   else          {max = first; min = first + step*(n-1);}

376:   is->min  = n > 0 ? min : PETSC_MAX_INT;
377:   is->max  = n > 0 ? max : PETSC_MIN_INT;
378:   is->data = (void*)sub;

380:   if ((!first && step == 1) || (first == max && step == -1 && !min)) is->isperm = PETSC_TRUE;
381:   else is->isperm = PETSC_FALSE;
382:   is->isidentity = PETSC_FALSE;
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 on MPI_Comm

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:   Concepts: IS^stride
409:   Concepts: index sets^stride
410:   Concepts: stride^index set

412: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
413: @*/
414: PetscErrorCode  ISCreateStride(MPI_Comm comm,PetscInt n,PetscInt first,PetscInt step,IS *is)
415: {

419:   ISCreate(comm,is);
420:   ISSetType(*is,ISSTRIDE);
421:   ISStrideSetStride(*is,n,first,step);
422:   return(0);
423: }

425: PETSC_EXTERN PetscErrorCode ISCreate_Stride(IS is)
426: {
428:   IS_Stride      *sub;

431:   PetscNewLog(is,&sub);
432:   is->data = (void *) sub;
433:   PetscMemcpy(is->ops,&myops,sizeof(myops));
434:   PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",ISStrideSetStride_Stride);
435:   return(0);
436: }