Actual source code: stride.c

      1: #define PETSCVEC_DLL
      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 PetscErrorCode PETSCVEC_DLLEXPORT VecInitializePackage(char *);

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

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

     32: PetscErrorCode ISDuplicate_Stride(IS is,IS *newIS)
     33: {
     35:   IS_Stride *sub = (IS_Stride*)is->data;

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

     44: PetscErrorCode ISInvertPermutation_Stride(IS is,PetscInt nlocal,IS *perm)
     45: {
     46:   IS_Stride *isstride = (IS_Stride*)is->data;

     50:   if (is->isidentity) {
     51:     ISCreateStride(PETSC_COMM_SELF,isstride->n,0,1,perm);
     52:   } else {
     53:     IS  tmp;
     54:     PetscInt *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:
     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: PetscErrorCode PETSCVEC_DLLEXPORT ISStrideGetInfo(IS is,PetscInt *first,PetscInt *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: }

    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: PetscErrorCode PETSCVEC_DLLEXPORT ISStride(IS is,PetscTruth *flag)
    127: {

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

    135:   return(0);
    136: }

    140: PetscErrorCode ISDestroy_Stride(IS is)
    141: {

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

    150: /*
    151:      Returns a legitimate index memory even if
    152:    the stride index set is empty.
    153: */
    156: PetscErrorCode ISGetIndices_Stride(IS in,PetscInt **idx)
    157: {
    158:   IS_Stride *sub = (IS_Stride*)in->data;
    160:   PetscInt       i;

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

    173: PetscErrorCode ISRestoreIndices_Stride(IS in,PetscInt **idx)
    174: {

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

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

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

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

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

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

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

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

    258: PetscErrorCode ISSorted_Stride(IS is,PetscTruth* flg)
    259: {
    260:   IS_Stride *sub = (IS_Stride*)is->data;

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

    268: static struct _ISOps myops = { ISGetSize_Stride,
    269:                                ISGetLocalSize_Stride,
    270:                                ISGetIndices_Stride,
    271:                                ISRestoreIndices_Stride,
    272:                                ISInvertPermutation_Stride,
    273:                                ISSort_Stride,
    274:                                ISSorted_Stride,
    275:                                ISDuplicate_Stride,
    276:                                ISDestroy_Stride,
    277:                                ISView_Stride,
    278:                                ISIdentity_Stride };

    282: /*@C
    283:    ISCreateStride - Creates a data structure for an index set
    284:    containing a list of evenly spaced integers.

    286:    Collective on MPI_Comm

    288:    Input Parameters:
    289: +  comm - the MPI communicator
    290: .  n - the length of the index set
    291: .  first - the first element of the index set
    292: -  step - the change to the next index

    294:    Output Parameter:
    295: .  is - the new index set

    297:    Notes:
    298:    When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
    299:    conceptually the same as MPI_Group operations. The IS are the
    300:    distributed sets of indices and thus certain operations on them are collective.

    302:    Level: beginner

    304:   Concepts: IS^stride
    305:   Concepts: index sets^stride
    306:   Concepts: stride^index set

    308: .seealso: ISCreateGeneral(), ISCreateBlock(), ISAllGather()
    309: @*/
    310: PetscErrorCode PETSCVEC_DLLEXPORT ISCreateStride(MPI_Comm comm,PetscInt n,PetscInt first,PetscInt step,IS *is)
    311: {
    313:   PetscInt        min,max;
    314:   IS         Nindex;
    315:   IS_Stride  *sub;
    316:   PetscTruth flg;

    320:   *is = PETSC_NULL;
    321:   if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Number of indices < 0");
    322: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
    323:   VecInitializePackage(PETSC_NULL);
    324: #endif

    326:   PetscHeaderCreate(Nindex,_p_IS,struct _ISOps,IS_COOKIE,IS_STRIDE,"IS",comm,ISDestroy,ISView);
    327:   PetscLogObjectMemory(Nindex,sizeof(IS_Stride) + sizeof(struct _p_IS));
    328:   PetscNew(IS_Stride,&sub);
    329:   sub->n         = n;
    330:   MPI_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,comm);
    331:   sub->first     = first;
    332:   sub->step      = step;
    333:   if (step > 0) {min = first; max = first + step*(n-1);}
    334:   else          {max = first; min = first + step*(n-1);}

    336:   Nindex->min     = min;
    337:   Nindex->max     = max;
    338:   Nindex->data    = (void*)sub;
    339:   PetscMemcpy(Nindex->ops,&myops,sizeof(myops));

    341:   if ((!first && step == 1) || (first == max && step == -1 && !min)) {
    342:     Nindex->isperm  = PETSC_TRUE;
    343:   } else {
    344:     Nindex->isperm  = PETSC_FALSE;
    345:   }
    346:   PetscOptionsHasName(PETSC_NULL,"-is_view",&flg);
    347:   if (flg) {
    348:     ISView(Nindex,PETSC_VIEWER_STDOUT_(Nindex->comm));
    349:   }
    350:   *is = Nindex;
    351:   return(0);
    352: }