Actual source code: stride.c

petsc-3.5.1 2014-08-06
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>             /*I   "petscis.h"   I*/
  7: #include <petscvec.h>
  8: #include <petscviewer.h>

 10: typedef struct {
 11:   PetscInt N,n,first,step;
 12: } IS_Stride;

 16: PetscErrorCode ISIdentity_Stride(IS is,PetscBool  *ident)
 17: {
 18:   IS_Stride *is_stride = (IS_Stride*)is->data;

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

 32: static PetscErrorCode ISCopy_Stride(IS is,IS isy)
 33: {
 34:   IS_Stride      *is_stride = (IS_Stride*)is->data,*isy_stride = (IS_Stride*)isy->data;

 38:   PetscMemcpy(isy_stride,is_stride,sizeof(IS_Stride));
 39:   return(0);
 40: }

 44: PetscErrorCode ISDuplicate_Stride(IS is,IS *newIS)
 45: {
 47:   IS_Stride      *sub = (IS_Stride*)is->data;

 50:   ISCreateStride(PetscObjectComm((PetscObject)is),sub->n,sub->first,sub->step,newIS);
 51:   return(0);
 52: }

 56: PetscErrorCode ISInvertPermutation_Stride(IS is,PetscInt nlocal,IS *perm)
 57: {
 58:   IS_Stride      *isstride = (IS_Stride*)is->data;

 62:   if (is->isidentity) {
 63:     ISCreateStride(PETSC_COMM_SELF,isstride->n,0,1,perm);
 64:   } else {
 65:     IS             tmp;
 66:     const PetscInt *indices,n = isstride->n;
 67:     ISGetIndices(is,&indices);
 68:     ISCreateGeneral(PetscObjectComm((PetscObject)is),n,indices,PETSC_COPY_VALUES,&tmp);
 69:     ISSetPermutation(tmp);
 70:     ISRestoreIndices(is,&indices);
 71:     ISInvertPermutation(tmp,nlocal,perm);
 72:     ISDestroy(&tmp);
 73:   }
 74:   return(0);
 75: }

 79: /*@
 80:    ISStrideGetInfo - Returns the first index in a stride index set and
 81:    the stride width.

 83:    Not Collective

 85:    Input Parameter:
 86: .  is - the index set

 88:    Output Parameters:
 89: .  first - the first index
 90: .  step - the stride width

 92:    Level: intermediate

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

 98:    Concepts: index sets^getting information
 99:    Concepts: IS^getting information

101: .seealso: ISCreateStride(), ISGetSize()
102: @*/
103: PetscErrorCode  ISStrideGetInfo(IS is,PetscInt *first,PetscInt *step)
104: {
105:   IS_Stride      *sub;
106:   PetscBool      flg;

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

116:   sub = (IS_Stride*)is->data;
117:   if (first) *first = sub->first;
118:   if (step)  *step  = sub->step;
119:   return(0);
120: }

124: PetscErrorCode ISDestroy_Stride(IS is)
125: {

129:   PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",0);
130:   PetscFree(is->data);
131:   return(0);
132: }

136: PetscErrorCode  ISToGeneral_Stride(IS inis)
137: {
139:   const PetscInt *idx;
140:   PetscInt       n;

143:   ISGetLocalSize(inis,&n);
144:   ISGetIndices(inis,&idx);
145:   ISSetType(inis,ISGENERAL);
146:   ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
147:   return(0);
148: }


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

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

174: PetscErrorCode ISRestoreIndices_Stride(IS in,const PetscInt *idx[])
175: {

179:   PetscFree(*(void**)idx);
180:   return(0);
181: }

185: PetscErrorCode ISGetSize_Stride(IS is,PetscInt *size)
186: {
187:   IS_Stride *sub = (IS_Stride*)is->data;

190:   *size = sub->N;
191:   return(0);
192: }

196: PetscErrorCode ISGetLocalSize_Stride(IS is,PetscInt *size)
197: {
198:   IS_Stride *sub = (IS_Stride*)is->data;

201:   *size = sub->n;
202:   return(0);
203: }

207: PetscErrorCode ISView_Stride(IS is,PetscViewer viewer)
208: {
209:   IS_Stride      *sub = (IS_Stride*)is->data;
210:   PetscInt       i,n = sub->n;
211:   PetscMPIInt    rank,size;
212:   PetscBool      iascii;

216:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
217:   if (iascii) {
218:     MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
219:     MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
220:     if (size == 1) {
221:       if (is->isperm) {
222:         PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");
223:       }
224:       PetscViewerASCIIPrintf(viewer,"Number of indices in (stride) set %D\n",n);
225:       for (i=0; i<n; i++) {
226:         PetscViewerASCIIPrintf(viewer,"%D %D\n",i,sub->first + i*sub->step);
227:       }
228:       PetscViewerFlush(viewer);
229:     } else {
230:       PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
231:       if (is->isperm) {
232:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutation\n",rank);
233:       }
234:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in (stride) set %D\n",rank,n);
235:       for (i=0; i<n; i++) {
236:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,sub->first + i*sub->step);
237:       }
238:       PetscViewerFlush(viewer);
239:       PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
240:     }
241:   }
242:   return(0);
243: }

247: PetscErrorCode ISSort_Stride(IS is)
248: {
249:   IS_Stride *sub = (IS_Stride*)is->data;

252:   if (sub->step >= 0) return(0);
253:   sub->first += (sub->n - 1)*sub->step;
254:   sub->step  *= -1;
255:   return(0);
256: }

260: PetscErrorCode ISSorted_Stride(IS is,PetscBool * flg)
261: {
262:   IS_Stride *sub = (IS_Stride*)is->data;

265:   if (sub->step >= 0) *flg = PETSC_TRUE;
266:   else *flg = PETSC_FALSE;
267:   return(0);
268: }

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,sub->n,sub->first,sub->step,newis);
279:   return(0);
280: }

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

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


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

352: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
353: @*/
354: PetscErrorCode  ISStrideSetStride(IS is,PetscInt n,PetscInt first,PetscInt step)
355: {

359:   if (n < 0) SETERRQ1(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Negative length %d not valid", n);
360:   PetscUseMethod(is,"ISStrideSetStride_C",(IS,PetscInt,PetscInt,PetscInt),(is,n,first,step));
361:   return(0);
362: }

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

373:   sub->n     = n;
374:   MPI_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)is));
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 : -1;
381:   is->max  = n > 0 ? max : -1;
382:   is->data = (void*)sub;

384:   if ((!first && step == 1) || (first == max && step == -1 && !min)) is->isperm = PETSC_TRUE;
385:   else is->isperm = PETSC_FALSE;
386:   return(0);
387: }

391: /*@
392:    ISCreateStride - Creates a data structure for an index set
393:    containing a list of evenly spaced integers.

395:    Collective on MPI_Comm

397:    Input Parameters:
398: +  comm - the MPI communicator
399: .  n - the length of the locally owned portion of the index set
400: .  first - the first element of the locally owned portion of the index set
401: -  step - the change to the next index

403:    Output Parameter:
404: .  is - the new index set

406:    Notes:
407:    When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
408:    conceptually the same as MPI_Group operations. The IS are the
409:    distributed sets of indices and thus certain operations on them are collective.

411:    Level: beginner

413:   Concepts: IS^stride
414:   Concepts: index sets^stride
415:   Concepts: stride^index set

417: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
418: @*/
419: PetscErrorCode  ISCreateStride(MPI_Comm comm,PetscInt n,PetscInt first,PetscInt step,IS *is)
420: {

424:   ISCreate(comm,is);
425:   ISSetType(*is,ISSTRIDE);
426:   ISStrideSetStride(*is,n,first,step);
427:   return(0);
428: }

432: PETSC_EXTERN PetscErrorCode ISCreate_Stride(IS is)
433: {
435:   IS_Stride      *sub;

438:   PetscMemcpy(is->ops,&myops,sizeof(myops));
439:   PetscNewLog(is,&sub);
440:   is->data = sub;
441:   PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",ISStrideSetStride_Stride);
442:   return(0);
443: }