Actual source code: aobasic.c

      1: #define PETSCDM_DLL

      3: /*
      4:     The most basic AO application ordering routines. These store the
      5:   entire orderings on each processor.
      6: */

     8:  #include src/dm/ao/aoimpl.h
     9:  #include petscsys.h

     11: typedef struct {
     12:   PetscInt N;
     13:   PetscInt *app,*petsc;  /* app[i] is the partner for the ith PETSc slot */
     14:                          /* petsc[j] is the partner for the jth app slot */
     15: } AO_Basic;

     17: /*
     18:        All processors have the same data so processor 1 prints it
     19: */
     22: PetscErrorCode AOView_Basic(AO ao,PetscViewer viewer)
     23: {
     25:   PetscMPIInt    rank;
     26:   PetscInt       i;
     27:   AO_Basic       *aodebug = (AO_Basic*)ao->data;
     28:   PetscTruth     iascii;

     31:   MPI_Comm_rank(ao->comm,&rank);
     32:   if (!rank){
     33:     PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
     34:     if (iascii) {
     35:       PetscViewerASCIIPrintf(viewer,"Number of elements in ordering %D\n",aodebug->N);
     36:       PetscViewerASCIIPrintf(viewer,  "PETSc->App  App->PETSc\n");
     37:       for (i=0; i<aodebug->N; i++) {
     38:         PetscViewerASCIIPrintf(viewer,"%3D  %3D    %3D  %3D\n",i,aodebug->app[i],i,aodebug->petsc[i]);
     39:       }
     40:     } else {
     41:       SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for AOData basic",((PetscObject)viewer)->type_name);
     42:     }
     43:   }
     44:   PetscViewerFlush(viewer);
     45:   return(0);
     46: }

     50: PetscErrorCode AODestroy_Basic(AO ao)
     51: {
     52:   AO_Basic       *aodebug = (AO_Basic*)ao->data;

     56:   PetscFree(aodebug->app);
     57:   PetscFree(ao->data);
     58:   return(0);
     59: }

     63: PetscErrorCode AOBasicGetIndices_Private(AO ao,PetscInt **app,PetscInt **petsc)
     64: {
     65:   AO_Basic *basic = (AO_Basic*)ao->data;

     68:   if (app)   *app   = basic->app;
     69:   if (petsc) *petsc = basic->petsc;
     70:   return(0);
     71: }

     75: PetscErrorCode AOPetscToApplication_Basic(AO ao,PetscInt n,PetscInt *ia)
     76: {
     77:   PetscInt i;
     78:   AO_Basic *aodebug = (AO_Basic*)ao->data;

     81:   for (i=0; i<n; i++) {
     82:     if (ia[i] >= 0) {ia[i] = aodebug->app[ia[i]];}
     83:   }
     84:   return(0);
     85: }

     89: PetscErrorCode AOApplicationToPetsc_Basic(AO ao,PetscInt n,PetscInt *ia)
     90: {
     91:   PetscInt i;
     92:   AO_Basic *aodebug = (AO_Basic*)ao->data;

     95:   for (i=0; i<n; i++) {
     96:     if (ia[i] >= 0) {ia[i] = aodebug->petsc[ia[i]];}
     97:   }
     98:   return(0);
     99: }

    103: PetscErrorCode AOPetscToApplicationPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
    104: {
    105:   AO_Basic       *aodebug = (AO_Basic *) ao->data;
    106:   PetscInt       *temp;
    107:   PetscInt       i, j;

    111:   PetscMalloc(aodebug->N*block * sizeof(PetscInt), &temp);
    112:   for(i = 0; i < aodebug->N; i++) {
    113:     for(j = 0; j < block; j++) temp[i*block+j] = array[aodebug->petsc[i]*block+j];
    114:   }
    115:   PetscMemcpy(array, temp, aodebug->N*block * sizeof(PetscInt));
    116:   PetscFree(temp);
    117:   return(0);
    118: }

    122: PetscErrorCode AOApplicationToPetscPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array)
    123: {
    124:   AO_Basic       *aodebug = (AO_Basic *) ao->data;
    125:   PetscInt       *temp;
    126:   PetscInt       i, j;

    130:   PetscMalloc(aodebug->N*block * sizeof(PetscInt), &temp);
    131:   for(i = 0; i < aodebug->N; i++) {
    132:     for(j = 0; j < block; j++) temp[i*block+j] = array[aodebug->app[i]*block+j];
    133:   }
    134:   PetscMemcpy(array, temp, aodebug->N*block * sizeof(PetscInt));
    135:   PetscFree(temp);
    136:   return(0);
    137: }

    141: PetscErrorCode AOPetscToApplicationPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
    142: {
    143:   AO_Basic       *aodebug = (AO_Basic *) ao->data;
    144:   PetscReal      *temp;
    145:   PetscInt       i, j;

    149:   PetscMalloc(aodebug->N*block * sizeof(PetscReal), &temp);
    150:   for(i = 0; i < aodebug->N; i++) {
    151:     for(j = 0; j < block; j++) temp[i*block+j] = array[aodebug->petsc[i]*block+j];
    152:   }
    153:   PetscMemcpy(array, temp, aodebug->N*block * sizeof(PetscReal));
    154:   PetscFree(temp);
    155:   return(0);
    156: }

    160: PetscErrorCode AOApplicationToPetscPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array)
    161: {
    162:   AO_Basic       *aodebug = (AO_Basic *) ao->data;
    163:   PetscReal      *temp;
    164:   PetscInt       i, j;

    168:   PetscMalloc(aodebug->N*block * sizeof(PetscReal), &temp);
    169:   for(i = 0; i < aodebug->N; i++) {
    170:     for(j = 0; j < block; j++) temp[i*block+j] = array[aodebug->app[i]*block+j];
    171:   }
    172:   PetscMemcpy(array, temp, aodebug->N*block * sizeof(PetscReal));
    173:   PetscFree(temp);
    174:   return(0);
    175: }

    177: static struct _AOOps AOops = {AOView_Basic,
    178:                               AODestroy_Basic,
    179:                               AOPetscToApplication_Basic,
    180:                               AOApplicationToPetsc_Basic,
    181:                               AOPetscToApplicationPermuteInt_Basic,
    182:                               AOApplicationToPetscPermuteInt_Basic,
    183:                               AOPetscToApplicationPermuteReal_Basic,
    184:                               AOApplicationToPetscPermuteReal_Basic};

    188: /*@C
    189:    AOCreateBasic - Creates a basic application ordering using two integer arrays.

    191:    Collective on MPI_Comm

    193:    Input Parameters:
    194: +  comm - MPI communicator that is to share AO
    195: .  napp - size of integer arrays
    196: .  myapp - integer array that defines an ordering
    197: -  mypetsc - integer array that defines another ordering (may be PETSC_NULL to
    198:              indicate the natural ordering)

    200:    Output Parameter:
    201: .  aoout - the new application ordering

    203:    Options Database Key:
    204: .   -ao_view - call AOView() at the conclusion of AOCreateBasic()

    206:    Level: beginner

    208: .keywords: AO, create

    210: .seealso: AOCreateBasicIS(), AODestroy()
    211: @*/
    212: PetscErrorCode PETSCDM_DLLEXPORT AOCreateBasic(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
    213: {
    214:   AO_Basic       *aobasic;
    215:   AO             ao;
    216:   PetscMPIInt    *lens,size,rank,nnapp,*disp;
    217:   PetscInt       *allpetsc,*allapp,ip,ia,N,i,*petsc,start;
    218:   PetscTruth     opt;

    223:   *aoout = 0;
    224: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
    225:   DMInitializePackage(PETSC_NULL);
    226: #endif

    228:   PetscHeaderCreate(ao, _p_AO, struct _AOOps, AO_COOKIE, AO_BASIC, "AO", comm, AODestroy, AOView);
    229:   PetscNew(AO_Basic, &aobasic);
    230:   PetscLogObjectMemory(ao, sizeof(struct _p_AO) + sizeof(AO_Basic));

    232:   PetscMemcpy(ao->ops, &AOops, sizeof(AOops));
    233:   ao->data = (void*) aobasic;

    235:   /* transmit all lengths to all processors */
    236:   MPI_Comm_size(comm, &size);
    237:   MPI_Comm_rank(comm, &rank);
    238:   PetscMalloc(2*size * sizeof(PetscMPIInt), &lens);
    239:   disp  = lens + size;
    240:   nnapp = napp;
    241:   MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm);
    242:   N    =  0;
    243:   for(i = 0; i < size; i++) {
    244:     disp[i] = N;
    245:     N += lens[i];
    246:   }
    247:   aobasic->N = N;

    249:   /*
    250:      If mypetsc is 0 then use "natural" numbering
    251:   */
    252:   if (napp && !mypetsc) {
    253:     start = disp[rank];
    254:     PetscMalloc((napp+1) * sizeof(PetscInt), &petsc);
    255:     for (i=0; i<napp; i++) {
    256:       petsc[i] = start + i;
    257:     }
    258:   } else {
    259:     petsc = (PetscInt*)mypetsc;
    260:   }

    262:   /* get all indices on all processors */
    263:   PetscMalloc(2*N * sizeof(PetscInt), &allpetsc);
    264:   allapp = allpetsc + N;
    265:   MPI_Allgatherv(petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
    266:   MPI_Allgatherv((void*)myapp, napp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm);
    267:   PetscFree(lens);

    269: #if defined(PETSC_USE_DEBUG)
    270:   {
    271:     PetscInt *sorted;
    272:     PetscMalloc((N+1)*sizeof(PetscInt),&sorted);

    274:     PetscMemcpy(sorted,allapp,N*sizeof(PetscInt));
    275:     PetscSortInt(N,sorted);
    276:     for (i=0; i<N; i++) {
    277:       if (sorted[i] != i) SETERRQ2(PETSC_ERR_ARG_WRONG,"PETSc ordering is missing %D has %D",i,sorted[i]);
    278:     }

    280:     PetscMemcpy(sorted,allapp,N*sizeof(PetscInt));
    281:     PetscSortInt(N,sorted);
    282:     for (i=0; i<N; i++) {
    283:       if (sorted[i] != i) SETERRQ2(PETSC_ERR_ARG_WRONG,"Application ordering is missing %D has %D",i,sorted[i]);
    284:     }

    286:     PetscFree(sorted);
    287:   }
    288: #endif

    290:   /* generate a list of application and PETSc node numbers */
    291:   PetscMalloc(2*N * sizeof(PetscInt), &aobasic->app);
    292:   PetscLogObjectMemory(ao,2*N*sizeof(PetscInt));
    293:   aobasic->petsc = aobasic->app + N;
    294:   PetscMemzero(aobasic->app, 2*N*sizeof(PetscInt));
    295:   for(i = 0; i < N; i++) {
    296:     ip = allpetsc[i];
    297:     ia = allapp[i];
    298:     /* check there are no duplicates */
    299:     if (aobasic->app[ip]) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"Duplicate in PETSc ordering at position %d. Already mapped to %d, not %d.", i, aobasic->app[ip]-1, ia);
    300:     aobasic->app[ip] = ia + 1;
    301:     if (aobasic->petsc[ia]) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"Duplicate in Application ordering at position %d. Already mapped to %d, not %d.", i, aobasic->petsc[ia]-1, ip);
    302:     aobasic->petsc[ia] = ip + 1;
    303:   }
    304:   if (!mypetsc) {
    305:     PetscFree(petsc);
    306:   }
    307:   PetscFree(allpetsc);
    308:   /* shift indices down by one */
    309:   for(i = 0; i < N; i++) {
    310:     aobasic->app[i]--;
    311:     aobasic->petsc[i]--;
    312:   }

    314:   PetscOptionsHasName(PETSC_NULL, "-ao_view", &opt);
    315:   if (opt) {
    316:     AOView(ao, PETSC_VIEWER_STDOUT_SELF);
    317:   }

    319:   *aoout = ao;
    320:   return(0);
    321: }

    325: /*@C
    326:    AOCreateBasicIS - Creates a basic application ordering using two index sets.

    328:    Collective on IS

    330:    Input Parameters:
    331: +  isapp - index set that defines an ordering
    332: -  ispetsc - index set that defines another ordering (may be PETSC_NULL to use the
    333:              natural ordering)

    335:    Output Parameter:
    336: .  aoout - the new application ordering

    338:    Options Database Key:
    339: -   -ao_view - call AOView() at the conclusion of AOCreateBasicIS()

    341:    Level: beginner

    343: .keywords: AO, create

    345: .seealso: AOCreateBasic(),  AODestroy()
    346: @*/
    347: PetscErrorCode PETSCDM_DLLEXPORT AOCreateBasicIS(IS isapp,IS ispetsc,AO *aoout)
    348: {
    350:   PetscInt       *mypetsc = 0,*myapp,napp,npetsc;
    351:   MPI_Comm       comm;

    354:   PetscObjectGetComm((PetscObject)isapp,&comm);
    355:   ISGetLocalSize(isapp,&napp);
    356:   if (ispetsc) {
    357:     ISGetLocalSize(ispetsc,&npetsc);
    358:     if (napp != npetsc) SETERRQ(PETSC_ERR_ARG_SIZ,"Local IS lengths must match");
    359:     ISGetIndices(ispetsc,&mypetsc);
    360:   }
    361:   ISGetIndices(isapp,&myapp);

    363:   AOCreateBasic(comm,napp,myapp,mypetsc,aoout);

    365:   ISRestoreIndices(isapp,&myapp);
    366:   if (ispetsc) {
    367:     ISRestoreIndices(ispetsc,&mypetsc);
    368:   }
    369:   return(0);
    370: }