Actual source code: general.c

      1: #define PETSCVEC_DLL
      2: /*
      3:      Provides the functions for index sets (IS) defined by a list of integers.
      4: */
     5:  #include src/vec/is/impls/general/general.h

      7: EXTERN PetscErrorCode PETSCVEC_DLLEXPORT VecInitializePackage(char *);

     11: PetscErrorCode ISDuplicate_General(IS is,IS *newIS)
     12: {
     14:   IS_General     *sub = (IS_General *)is->data;

     17:   ISCreateGeneral(is->comm,sub->n,sub->idx,newIS);
     18:   return(0);
     19: }

     23: PetscErrorCode ISDestroy_General(IS is)
     24: {
     25:   IS_General     *is_general = (IS_General*)is->data;

     29:   if (is_general->allocated) {
     30:     PetscFree(is_general->idx);
     31:   }
     32:   PetscFree(is_general);
     33:   PetscHeaderDestroy(is);
     34:   return(0);
     35: }

     39: PetscErrorCode ISIdentity_General(IS is,PetscTruth *ident)
     40: {
     41:   IS_General *is_general = (IS_General*)is->data;
     42:   PetscInt   i,n = is_general->n,*idx = is_general->idx;

     45:   is->isidentity = PETSC_TRUE;
     46:   *ident         = PETSC_TRUE;
     47:   for (i=0; i<n; i++) {
     48:     if (idx[i] != i) {
     49:       is->isidentity = PETSC_FALSE;
     50:       *ident         = PETSC_FALSE;
     51:       break;
     52:     }
     53:   }
     54:   return(0);
     55: }

     59: PetscErrorCode ISGetIndices_General(IS in,PetscInt **idx)
     60: {
     61:   IS_General *sub = (IS_General*)in->data;

     64:   *idx = sub->idx;
     65:   return(0);
     66: }

     70: PetscErrorCode ISRestoreIndices_General(IS in,PetscInt **idx)
     71: {
     72:   IS_General *sub = (IS_General*)in->data;

     75:   if (*idx != sub->idx) {
     76:     SETERRQ(PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
     77:   }
     78:   return(0);
     79: }

     83: PetscErrorCode ISGetSize_General(IS is,PetscInt *size)
     84: {
     85:   IS_General *sub = (IS_General *)is->data;

     88:   *size = sub->N;
     89:   return(0);
     90: }

     94: PetscErrorCode ISGetLocalSize_General(IS is,PetscInt *size)
     95: {
     96:   IS_General *sub = (IS_General *)is->data;

     99:   *size = sub->n;
    100:   return(0);
    101: }

    105: PetscErrorCode ISInvertPermutation_General(IS is,PetscInt nlocal,IS *isout)
    106: {
    107:   IS_General     *sub = (IS_General *)is->data;
    108:   PetscInt       i,*ii,n = sub->n,*idx = sub->idx,nstart;
    109:   PetscMPIInt    size;
    110:   IS             istmp,nistmp;

    114:   MPI_Comm_size(is->comm,&size);
    115:   if (size == 1) {
    116:     PetscMalloc(n*sizeof(PetscInt),&ii);
    117:     for (i=0; i<n; i++) {
    118:       ii[idx[i]] = i;
    119:     }
    120:     ISCreateGeneral(PETSC_COMM_SELF,n,ii,isout);
    121:     ISSetPermutation(*isout);
    122:     PetscFree(ii);
    123:   } else {
    124:     /* crude, nonscalable get entire IS on each processor */
    125:     if (nlocal == PETSC_DECIDE) SETERRQ(PETSC_ERR_SUP,"Do not yet support nlocal of PETSC_DECIDE");
    126:     ISAllGather(is,&istmp);
    127:     ISSetPermutation(istmp);
    128:     ISInvertPermutation(istmp,PETSC_DECIDE,&nistmp);
    129:     ISDestroy(istmp);
    130:     /* get the part we need */
    131:     MPI_Scan(&nlocal,&nstart,1,MPIU_INT,MPI_SUM,is->comm);
    132: #if defined(PETSC_USE_DEBUG)
    133:     {
    134:       PetscMPIInt rank;
    135:       MPI_Comm_rank(is->comm,&rank);
    136:       if (rank == size-1) {
    137:         if (nstart != sub->N) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Sum of nlocal lengths %d != total IS length %d",nstart,sub->N);
    138:       }
    139:     }
    140: #endif
    141:     nstart -= nlocal;
    142:     ISGetIndices(nistmp,&idx);
    143:     ISCreateGeneral(is->comm,nlocal,idx+nstart,isout);
    144:     ISRestoreIndices(nistmp,&idx);
    145:     ISDestroy(nistmp);
    146:   }
    147:   return(0);
    148: }

    152: PetscErrorCode ISView_General(IS is,PetscViewer viewer)
    153: {
    154:   IS_General     *sub = (IS_General *)is->data;
    156:   PetscInt       i,n = sub->n,*idx = sub->idx;
    157:   PetscTruth     iascii;

    160:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    161:   if (iascii) {
    162:     MPI_Comm    comm;
    163:     PetscMPIInt rank,size;

    165:     PetscObjectGetComm((PetscObject)viewer,&comm);
    166:     MPI_Comm_rank(comm,&rank);
    167:     MPI_Comm_size(comm,&size);

    169:     if (size > 1) {
    170:       if (is->isperm) {
    171:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutation\n",rank);
    172:       }
    173:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in set %D\n",rank,n);
    174:       for (i=0; i<n; i++) {
    175:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i,idx[i]);
    176:       }
    177:     } else {
    178:       if (is->isperm) {
    179:         PetscViewerASCIISynchronizedPrintf(viewer,"Index set is permutation\n");
    180:       }
    181:       PetscViewerASCIISynchronizedPrintf(viewer,"Number of indices in set %D\n",n);
    182:       for (i=0; i<n; i++) {
    183:         PetscViewerASCIISynchronizedPrintf(viewer,"%D %D\n",i,idx[i]);
    184:       }
    185:     }
    186:     PetscViewerFlush(viewer);
    187:   } else {
    188:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for this object",((PetscObject)viewer)->type_name);
    189:   }
    190:   return(0);
    191: }

    195: PetscErrorCode ISSort_General(IS is)
    196: {
    197:   IS_General     *sub = (IS_General *)is->data;

    201:   if (sub->sorted) return(0);
    202:   PetscSortInt(sub->n,sub->idx);
    203:   sub->sorted = PETSC_TRUE;
    204:   return(0);
    205: }

    209: PetscErrorCode ISSorted_General(IS is,PetscTruth *flg)
    210: {
    211:   IS_General *sub = (IS_General *)is->data;

    214:   *flg = sub->sorted;
    215:   return(0);
    216: }

    218: static struct _ISOps myops = { ISGetSize_General,
    219:                                ISGetLocalSize_General,
    220:                                ISGetIndices_General,
    221:                                ISRestoreIndices_General,
    222:                                ISInvertPermutation_General,
    223:                                ISSort_General,
    224:                                ISSorted_General,
    225:                                ISDuplicate_General,
    226:                                ISDestroy_General,
    227:                                ISView_General,
    228:                                ISIdentity_General };

    232: PetscErrorCode ISCreateGeneral_Private(MPI_Comm comm,IS *is)
    233: {
    235:   IS             Nindex = *is;
    236:   IS_General     *sub = (IS_General*)Nindex->data;
    237:   PetscInt       n = sub->n,i,min,max;
    238:   const PetscInt *idx = sub->idx;
    239:   PetscTruth     sorted = PETSC_TRUE;
    240:   PetscTruth     flg;

    244:   if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
    246:   *is = PETSC_NULL;
    247: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
    248:   VecInitializePackage(PETSC_NULL);
    249: #endif

    251:   MPI_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,comm);
    252:   for (i=1; i<n; i++) {
    253:     if (idx[i] < idx[i-1]) {sorted = PETSC_FALSE; break;}
    254:   }
    255:   if (n) {min = max = idx[0];} else {min = max = 0;}
    256:   for (i=1; i<n; i++) {
    257:     if (idx[i] < min) min = idx[i];
    258:     if (idx[i] > max) max = idx[i];
    259:   }
    260:   sub->sorted     = sorted;
    261:   Nindex->min     = min;
    262:   Nindex->max     = max;
    263:   PetscMemcpy(Nindex->ops,&myops,sizeof(myops));
    264:   Nindex->isperm     = PETSC_FALSE;
    265:   Nindex->isidentity = PETSC_FALSE;
    266:   PetscOptionsHasName(PETSC_NULL,"-is_view",&flg);
    267:   if (flg) {
    268:     ISView(Nindex,PETSC_VIEWER_STDOUT_(Nindex->comm));
    269:   }
    270:   *is = Nindex;
    271:   return(0);
    272: }

    274: /*@C
    275:    ISCreateGeneral - Creates a data structure for an index set
    276:    containing a list of integers.

    278:    Collective on MPI_Comm

    280:    Input Parameters:
    281: +  comm - the MPI communicator
    282: .  n - the length of the index set
    283: -  idx - the list of integers

    285:    Output Parameter:
    286: .  is - the new index set

    288:    Notes:
    289:    The index array is copied to internally allocated storage. After the call,
    290:    the user can free the index array.

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

    297:    Level: beginner

    299:   Concepts: index sets^creating
    300:   Concepts: IS^creating

    302: .seealso: ISCreateGeneralWithArray(), ISCreateStride(), ISCreateBlock(), ISAllGather()
    303: @*/
    304: PetscErrorCode PETSCVEC_DLLEXPORT ISCreateGeneral(MPI_Comm comm,PetscInt n,const PetscInt idx[],IS *is)
    305: {
    307:   IS             Nindex;
    308:   IS_General     *sub;

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

    319:   PetscHeaderCreate(Nindex,_p_IS,struct _ISOps,IS_COOKIE,IS_GENERAL,"IS",comm,ISDestroy,ISView);
    320:   PetscNew(IS_General,&sub);
    321:   PetscLogObjectMemory(Nindex,sizeof(IS_General)+n*sizeof(PetscInt)+sizeof(struct _p_IS));
    322:   PetscMalloc(n*sizeof(PetscInt),&sub->idx);
    323:   PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));
    324:   sub->n         = n;
    325:   sub->allocated = PETSC_TRUE;
    326:   Nindex->data   = (void*)sub;

    328:   *is = Nindex;
    329:   ISCreateGeneral_Private(comm,is);

    331:   return(0);
    332: }

    334: /*@C
    335:    ISCreateGeneralWithArray - Creates a data structure for an index set
    336:    containing a list of integers.

    338:    Collective on MPI_Comm

    340:    Input Parameters:
    341: +  comm - the MPI communicator
    342: .  n - the length of the index set
    343: -  idx - the list of integers

    345:    Output Parameter:
    346: .  is - the new index set

    348:    Notes:
    349:    Unlike with ISCreateGeneral, the indices are not copied to internally
    350:    allocated storage. The user array is not freed by ISDestroy.

    352:    When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
    353:    conceptually the same as MPI_Group operations. The IS are then
    354:    distributed sets of indices and thus certain operations on them are collective.

    356:    Level: beginner

    358:   Concepts: index sets^creating
    359:   Concepts: IS^creating

    361: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
    362: @*/
    363: PetscErrorCode PETSCVEC_DLLEXPORT ISCreateGeneralWithArray(MPI_Comm comm,PetscInt n,PetscInt idx[],IS *is)
    364: {
    366:   IS             Nindex;
    367:   IS_General     *sub;

    371:   if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
    373:   *is = PETSC_NULL;
    374: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
    375:   VecInitializePackage(PETSC_NULL);
    376: #endif

    378:   PetscHeaderCreate(Nindex,_p_IS,struct _ISOps,IS_COOKIE,IS_GENERAL,"IS",comm,ISDestroy,ISView);
    379:   PetscNew(IS_General,&sub);
    380:   PetscLogObjectMemory(Nindex,sizeof(IS_General)+n*sizeof(PetscInt)+sizeof(struct _p_IS));
    381:   sub->idx       = idx;
    382:   sub->n         = n;
    383:   sub->allocated = PETSC_FALSE;
    384:   Nindex->data   = (void*)sub;

    386:   *is = Nindex;
    387:   ISCreateGeneral_Private(comm,is);

    389:   return(0);
    390: }