Actual source code: stride.c

petsc-master 2019-10-12
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: 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),is->map->n,sub->first,sub->step,newIS);
 44:   return(0);
 45: }

 47: PetscErrorCode ISInvertPermutation_Stride(IS is,PetscInt nlocal,IS *perm)
 48: {

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

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

 71:    Not Collective

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

 76:    Output Parameters:
 77: +  first - the first index
 78: -  step - the stride width

 80:    Level: intermediate

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


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

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

102:   sub = (IS_Stride*)is->data;
103:   if (first) *first = sub->first;
104:   if (step)  *step  = sub->step;
105:   return(0);
106: }

108: PetscErrorCode ISDestroy_Stride(IS is)
109: {

113:   PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",NULL);
114:   PetscFree(is->data);
115:   return(0);
116: }

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

125:   ISGetLocalSize(inis,&n);
126:   ISGetIndices(inis,&idx);
127:   ISSetType(inis,ISGENERAL);
128:   ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
129:   return(0);
130: }

132: PetscErrorCode ISLocate_Stride(IS is,PetscInt key,PetscInt *location)
133: {
134:   IS_Stride      *sub = (IS_Stride*)is->data;
135:   PetscInt       rem, step;

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

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

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

167: PetscErrorCode ISRestoreIndices_Stride(IS in,const PetscInt *idx[])
168: {

172:   PetscFree(*(void**)idx);
173:   return(0);
174: }

176: PetscErrorCode ISView_Stride(IS is,PetscViewer viewer)
177: {
178:   IS_Stride         *sub = (IS_Stride*)is->data;
179:   PetscInt          i,n = is->map->n;
180:   PetscMPIInt       rank,size;
181:   PetscBool         iascii;
182:   PetscViewerFormat fmt;
183:   PetscErrorCode    ierr;

186:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
187:   if (iascii) {
188:     PetscBool matl;

190:     MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
191:     MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
192:     PetscViewerGetFormat(viewer,&fmt);
193:     matl = (PetscBool)(fmt == PETSC_VIEWER_ASCII_MATLAB);
194:     if (size == 1) {
195:       if (matl) {
196:         const char* name;

198:         PetscObjectGetName((PetscObject)is,&name);
199:         PetscViewerASCIIPrintf(viewer,"%s = [%D : %D : %D];\n",name,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
200:       } else {
201:         if (is->isperm) {
202:           PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");
203:         }
204:         PetscViewerASCIIPrintf(viewer,"Number of indices in (stride) set %D\n",n);
205:         for (i=0; i<n; i++) {
206:           PetscViewerASCIIPrintf(viewer,"%D %D\n",i,sub->first + i*sub->step);
207:         }
208:       }
209:       PetscViewerFlush(viewer);
210:     } else {
211:       PetscViewerASCIIPushSynchronized(viewer);
212:       if (matl) {
213:         const char* name;

215:         PetscObjectGetName((PetscObject)is,&name);
216:         PetscViewerASCIISynchronizedPrintf(viewer,"%s_%d = [%D : %D : %D];\n",name,rank,sub->first+1,sub->step,sub->first + sub->step*(n-1)+1);
217:       } else {
218:         if (is->isperm) {
219:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutation\n",rank);
220:         }
221:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in (stride) set %D\n",rank,n);
222:         for (i=0; i<n; i++) {
223:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,sub->first + i*sub->step);
224:         }
225:       }
226:       PetscViewerFlush(viewer);
227:       PetscViewerASCIIPopSynchronized(viewer);
228:     }
229:   }
230:   return(0);
231: }

233: PetscErrorCode ISSort_Stride(IS is)
234: {
235:   IS_Stride *sub = (IS_Stride*)is->data;

238:   if (sub->step >= 0) return(0);
239:   sub->first += (is->map->n - 1)*sub->step;
240:   sub->step  *= -1;
241:   return(0);
242: }

244: PetscErrorCode ISSorted_Stride(IS is,PetscBool * flg)
245: {
246:   IS_Stride *sub = (IS_Stride*)is->data;

249:   if (sub->step >= 0) *flg = PETSC_TRUE;
250:   else *flg = PETSC_FALSE;
251:   return(0);
252: }

254: static PetscErrorCode ISOnComm_Stride(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
255: {
257:   IS_Stride      *sub = (IS_Stride*)is->data;

260:   ISCreateStride(comm,is->map->n,sub->first,sub->step,newis);
261:   return(0);
262: }

264: static PetscErrorCode ISSetBlockSize_Stride(IS is,PetscInt bs)
265: {
266:   IS_Stride     *sub = (IS_Stride*)is->data;

270:   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);
271:   PetscLayoutSetBlockSize(is->map, bs);
272:   return(0);
273: }

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

280:   if (sub->step == 1 && sub->first >= gstart && sub->first+is->map->n <= gend) {
281:     *start  = sub->first - gstart;
282:     *contig = PETSC_TRUE;
283:   } else {
284:     *start  = -1;
285:     *contig = PETSC_FALSE;
286:   }
287:   return(0);
288: }


291: static struct _ISOps myops = { ISGetIndices_Stride,
292:                                ISRestoreIndices_Stride,
293:                                ISInvertPermutation_Stride,
294:                                ISSort_Stride,
295:                                ISSort_Stride,
296:                                ISSorted_Stride,
297:                                ISDuplicate_Stride,
298:                                ISDestroy_Stride,
299:                                ISView_Stride,
300:                                ISLoad_Default,
301:                                ISIdentity_Stride,
302:                                ISCopy_Stride,
303:                                ISToGeneral_Stride,
304:                                ISOnComm_Stride,
305:                                ISSetBlockSize_Stride,
306:                                ISContiguousLocal_Stride,
307:                                ISLocate_Stride};


310: /*@
311:    ISStrideSetStride - Sets the stride information for a stride index set.

313:    Collective on IS

315:    Input Parameters:
316: +  is - the index set
317: .  n - the length of the locally owned portion of the index set
318: .  first - the first element of the locally owned portion of the index set
319: -  step - the change to the next index

321:    Level: beginner

323: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
324: @*/
325: PetscErrorCode  ISStrideSetStride(IS is,PetscInt n,PetscInt first,PetscInt step)
326: {

330:   if (n < 0) SETERRQ1(PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Negative length %d not valid", n);
331:   PetscUseMethod(is,"ISStrideSetStride_C",(IS,PetscInt,PetscInt,PetscInt),(is,n,first,step));
332:   return(0);
333: }

335: PetscErrorCode  ISStrideSetStride_Stride(IS is,PetscInt n,PetscInt first,PetscInt step)
336: {
338:   PetscInt       min,max;
339:   IS_Stride      *sub = (IS_Stride*)is->data;
340:   PetscLayout    map;

343:   PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is),n,is->map->N,is->map->bs,&map);
344:   PetscLayoutDestroy(&is->map);
345:   is->map = map;

347:   sub->first = first;
348:   sub->step  = step;
349:   if (step > 0) {min = first; max = first + step*(n-1);}
350:   else          {max = first; min = first + step*(n-1);}

352:   is->min  = n > 0 ? min : PETSC_MAX_INT;
353:   is->max  = n > 0 ? max : PETSC_MIN_INT;
354:   is->data = (void*)sub;

356:   if ((!first && step == 1) || (first == max && step == -1 && !min)) is->isperm = PETSC_TRUE;
357:   else is->isperm = PETSC_FALSE;
358:   is->isidentity = PETSC_FALSE;
359:   return(0);
360: }

362: /*@
363:    ISCreateStride - Creates a data structure for an index set
364:    containing a list of evenly spaced integers.

366:    Collective

368:    Input Parameters:
369: +  comm - the MPI communicator
370: .  n - the length of the locally owned portion of the index set
371: .  first - the first element of the locally owned portion of the index set
372: -  step - the change to the next index

374:    Output Parameter:
375: .  is - the new index set

377:    Notes:
378:    When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
379:    conceptually the same as MPI_Group operations. The IS are the
380:    distributed sets of indices and thus certain operations on them are collective.

382:    Level: beginner

384: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
385: @*/
386: PetscErrorCode  ISCreateStride(MPI_Comm comm,PetscInt n,PetscInt first,PetscInt step,IS *is)
387: {

391:   ISCreate(comm,is);
392:   ISSetType(*is,ISSTRIDE);
393:   ISStrideSetStride(*is,n,first,step);
394:   return(0);
395: }

397: PETSC_EXTERN PetscErrorCode ISCreate_Stride(IS is)
398: {
400:   IS_Stride      *sub;

403:   PetscNewLog(is,&sub);
404:   is->data = (void *) sub;
405:   PetscMemcpy(is->ops,&myops,sizeof(myops));
406:   PetscObjectComposeFunction((PetscObject)is,"ISStrideSetStride_C",ISStrideSetStride_Stride);
407:   return(0);
408: }