Actual source code: stride.c

  1: /*$Id: stride.c,v 1.104 2001/06/21 21:15:53 bsmith Exp $*/
  2: /*
  3:        Index sets of evenly space integers, defined by a 
  4:     start, stride and length.
  5: */
 6:  #include src/vec/is/isimpl.h

  8: EXTERN int VecInitializePackage(char *);

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

 14: #undef __FUNCT__  
 16: int ISIdentity_Stride(IS is,PetscTruth *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: }

 30: #undef __FUNCT__  
 32: int ISDuplicate_Stride(IS is,IS *newIS)
 33: {
 34:   int       ierr;
 35:   IS_Stride *sub = (IS_Stride*)is->data;

 38:   ISCreateStride(is->comm,sub->n,sub->first,sub->step,newIS);
 39:   return(0);
 40: }

 42: #undef __FUNCT__  
 44: int ISInvertPermutation_Stride(IS is,int nlocal,IS *perm)
 45: {
 46:   IS_Stride *isstride = (IS_Stride*)is->data;
 47:   int       ierr;

 50:   if (is->isidentity) {
 51:     ISCreateStride(PETSC_COMM_SELF,isstride->n,0,1,perm);
 52:   } else {
 53:     IS  tmp;
 54:     int *indices,n = isstride->n;
 55:     ISGetIndices(is,&indices);
 56:     ISCreateGeneral(is->comm,n,indices,&tmp);
 57:     ISRestoreIndices(is,&indices);
 58:     ISInvertPermutation(tmp,nlocal,perm);
 59:     ISDestroy(tmp);
 60:   }
 61:   return(0);
 62: }
 63: 
 64: #undef __FUNCT__  
 66: /*@
 67:    ISStrideGetInfo - Returns the first index in a stride index set and 
 68:    the stride width.

 70:    Not Collective

 72:    Input Parameter:
 73: .  is - the index set

 75:    Output Parameters:
 76: .  first - the first index
 77: .  step - the stride width

 79:    Level: intermediate

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

 85:    Concepts: index sets^getting information
 86:    Concepts: IS^getting information

 88: .seealso: ISCreateStride(), ISGetSize()
 89: @*/
 90: int ISStrideGetInfo(IS is,int *first,int *step)
 91: {
 92:   IS_Stride *sub;


 99:   sub = (IS_Stride*)is->data;
100:   if (is->type != IS_STRIDE) return(0);
101:   if (first) *first = sub->first;
102:   if (step)  *step  = sub->step;
103:   return(0);
104: }

106: #undef __FUNCT__  
108: /*@C
109:    ISStride - Determines if an IS is based on a stride.

111:    Not Collective

113:    Input Parameter:
114: .  is - the index set

116:    Output Parameters:
117: .  flag - either PETSC_TRUE or PETSC_FALSE

119:    Level: intermediate

121:    Concepts: index sets^is it stride
122:    Concepts: IS^is it stride

124: .seealso: ISCreateStride(), ISGetSize()
125: @*/
126: int ISStride(IS is,PetscTruth *flag)
127: {

132:   if (is->type != IS_STRIDE) *flag = PETSC_FALSE;
133:   else                       *flag = PETSC_TRUE;

135:   return(0);
136: }

138: #undef __FUNCT__  
140: int ISDestroy_Stride(IS is)
141: {

145:   PetscFree(is->data);
146:   PetscLogObjectDestroy(is);
147:   PetscHeaderDestroy(is); return(0);
148: }

150: /*
151:      Returns a legitimate index memory even if 
152:    the stride index set is empty.
153: */
154: #undef __FUNCT__  
156: int ISGetIndices_Stride(IS in,int **idx)
157: {
158:   IS_Stride *sub = (IS_Stride*)in->data;
159:   int       i,ierr;

162:   ierr      = PetscMalloc((sub->n+1)*sizeof(int),idx);
163:   (*idx)[0] = sub->first;
164:   for (i=1; i<sub->n; i++) (*idx)[i] = (*idx)[i-1] + sub->step;
165:   return(0);
166: }

168: #undef __FUNCT__  
170: int ISRestoreIndices_Stride(IS in,int **idx)
171: {

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

179: #undef __FUNCT__  
181: int ISGetSize_Stride(IS is,int *size)
182: {
183:   IS_Stride *sub = (IS_Stride *)is->data;

186:   *size = sub->N;
187:   return(0);
188: }

190: #undef __FUNCT__  
192: int ISGetLocalSize_Stride(IS is,int *size)
193: {
194:   IS_Stride *sub = (IS_Stride *)is->data;

197:   *size = sub->n;
198:   return(0);
199: }

201: #undef __FUNCT__  
203: int ISView_Stride(IS is,PetscViewer viewer)
204: {
205:   IS_Stride   *sub = (IS_Stride *)is->data;
206:   int         i,n = sub->n,ierr,rank,size;
207:   PetscTruth  isascii;

210:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
211:   if (isascii) {
212:     MPI_Comm_rank(is->comm,&rank);
213:     MPI_Comm_size(is->comm,&size);
214:     if (size == 1) {
215:       if (is->isperm) {
216:         PetscViewerASCIISynchronizedPrintf(viewer,"Index set is permutationn");
217:       }
218:       PetscViewerASCIISynchronizedPrintf(viewer,"Number of indices in (stride) set %dn",n);
219:       for (i=0; i<n; i++) {
220:         PetscViewerASCIISynchronizedPrintf(viewer,"%d %dn",i,sub->first + i*sub->step);
221:       }
222:     } else {
223:       if (is->isperm) {
224:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutationn",rank);
225:       }
226:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in (stride) set %dn",rank,n);
227:       for (i=0; i<n; i++) {
228:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %d %dn",rank,i,sub->first + i*sub->step);
229:       }
230:     }
231:     PetscViewerFlush(viewer);
232:   } else {
233:     SETERRQ1(1,"Viewer type %s not supported for this object",((PetscObject)viewer)->type_name);
234:   }
235:   return(0);
236: }
237: 
238: #undef __FUNCT__  
240: int ISSort_Stride(IS is)
241: {
242:   IS_Stride *sub = (IS_Stride*)is->data;

245:   if (sub->step >= 0) return(0);
246:   sub->first += (sub->n - 1)*sub->step;
247:   sub->step *= -1;
248:   return(0);
249: }

251: #undef __FUNCT__  
253: int ISSorted_Stride(IS is,PetscTruth* flg)
254: {
255:   IS_Stride *sub = (IS_Stride*)is->data;

258:   if (sub->step >= 0) *flg = PETSC_TRUE;
259:   else *flg = PETSC_FALSE;
260:   return(0);
261: }

263: static struct _ISOps myops = { ISGetSize_Stride,
264:                                ISGetLocalSize_Stride,
265:                                ISGetIndices_Stride,
266:                                ISRestoreIndices_Stride,
267:                                ISInvertPermutation_Stride,
268:                                ISSort_Stride,
269:                                ISSorted_Stride,
270:                                ISDuplicate_Stride,
271:                                ISDestroy_Stride,
272:                                ISView_Stride,
273:                                ISIdentity_Stride };

275: #undef __FUNCT__  
277: /*@C
278:    ISCreateStride - Creates a data structure for an index set 
279:    containing a list of evenly spaced integers.

281:    Collective on MPI_Comm

283:    Input Parameters:
284: +  comm - the MPI communicator
285: .  n - the length of the index set
286: .  first - the first element of the index set
287: -  step - the change to the next index

289:    Output Parameter:
290: .  is - the new index set

292:    Notes: 
293:    When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
294:    conceptually the same as MPI_Group operations. The IS are the 
295:    distributed sets of indices and thus certain operations on them are collective. 

297:    Level: beginner

299:   Concepts: IS^stride
300:   Concepts: index sets^stride
301:   Concepts: stride^index set

303: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
304: @*/
305: int ISCreateStride(MPI_Comm comm,int n,int first,int step,IS *is)
306: {
307:   int        min,max,ierr;
308:   IS         Nindex;
309:   IS_Stride  *sub;
310:   PetscTruth flg;

314:   *is = PETSC_NULL;
315:   if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Number of indices < 0");
316: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
317:   VecInitializePackage(PETSC_NULL);
318: #endif

320:   PetscHeaderCreate(Nindex,_p_IS,struct _ISOps,IS_COOKIE,IS_STRIDE,"IS",comm,ISDestroy,ISView);
321:   PetscLogObjectCreate(Nindex);
322:   PetscLogObjectMemory(Nindex,sizeof(IS_Stride) + sizeof(struct _p_IS));
323:   ierr           = PetscNew(IS_Stride,&sub);
324:   sub->n         = n;
325:   MPI_Allreduce(&n,&sub->N,1,MPI_INT,MPI_SUM,comm);
326:   sub->first     = first;
327:   sub->step      = step;
328:   if (step > 0) {min = first; max = first + step*(n-1);}
329:   else          {max = first; min = first + step*(n-1);}

331:   Nindex->min     = min;
332:   Nindex->max     = max;
333:   Nindex->data    = (void*)sub;
334:   PetscMemcpy(Nindex->ops,&myops,sizeof(myops));

336:   if ((first == 0 && step == 1) || (first == max && step == -1 && min == 0)) {
337:     Nindex->isperm  = PETSC_TRUE;
338:   } else {
339:     Nindex->isperm  = PETSC_FALSE;
340:   }
341:   PetscOptionsHasName(PETSC_NULL,"-is_view",&flg);
342:   if (flg) {
343:     ISView(Nindex,PETSC_VIEWER_STDOUT_(Nindex->comm));
344:   }
345:   *is = Nindex;
346:   return(0);
347: }