Actual source code: stride.c

petsc-3.11.4 2019-09-28
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 <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: }