Actual source code: aomapping.c

      1: #define PETSCDM_DLL

      3: /*
      4:   These AO application ordering routines do not require that the input
      5:   be a permutation, but merely a 1-1 mapping. This implementation still
      6:   keeps the entire ordering on each processor.
      7: */

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

     12: typedef struct {
     13:   PetscInt N;
     14:   PetscInt *app;       /* app[i] is the partner for petsc[appPerm[i]] */
     15:   PetscInt *appPerm;
     16:   PetscInt *petsc;     /* petsc[j] is the partner for app[petscPerm[j]] */
     17:   PetscInt *petscPerm;
     18: } AO_Mapping;

     22: PetscErrorCode AODestroy_Mapping(AO ao)
     23: {
     24:   AO_Mapping     *aomap = (AO_Mapping *) ao->data;

     28:   PetscFree(aomap->app);
     29:   PetscFree(ao->data);
     30:   return(0);
     31: }

     35: PetscErrorCode AOView_Mapping(AO ao, PetscViewer viewer)
     36: {
     37:   AO_Mapping     *aomap = (AO_Mapping *) ao->data;
     38:   PetscMPIInt    rank;
     39:   PetscInt       i;
     40:   PetscTruth     iascii;

     44:   MPI_Comm_rank(ao->comm, &rank);
     45:   if (rank) return(0);

     47:   if (!viewer) {
     48:     viewer = PETSC_VIEWER_STDOUT_SELF;
     49:   }

     51:   PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_ASCII, &iascii);
     52:   if (iascii) {
     53:     PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %D\n", aomap->N);
     54:     PetscViewerASCIIPrintf(viewer, "   App.   PETSc\n");
     55:     for(i = 0; i < aomap->N; i++) {
     56:       PetscViewerASCIIPrintf(viewer, "%D   %D    %D\n", i, aomap->app[i], aomap->petsc[aomap->appPerm[i]]);
     57:     }
     58:   }
     59:   return(0);
     60: }

     64: PetscErrorCode AOPetscToApplication_Mapping(AO ao, PetscInt n, PetscInt *ia)
     65: {
     66:   AO_Mapping *aomap = (AO_Mapping *) ao->data;
     67:   PetscInt   *app   = aomap->app;
     68:   PetscInt   *petsc = aomap->petsc;
     69:   PetscInt   *perm  = aomap->petscPerm;
     70:   PetscInt   N     = aomap->N;
     71:   PetscInt   low, high, mid=0;
     72:   PetscInt   idex;
     73:   PetscInt   i;

     75:   /* It would be possible to use a single bisection search, which
     76:      recursively divided the indices to be converted, and searched
     77:      partitions which contained an index. This would result in
     78:      better running times if indices are clustered.
     79:   */
     81:   for(i = 0; i < n; i++) {
     82:     idex = ia[i];
     83:     if (idex < 0) continue;
     84:     /* Use bisection since the array is sorted */
     85:     low  = 0;
     86:     high = N - 1;
     87:     while (low <= high) {
     88:       mid = (low + high)/2;
     89:       if (idex == petsc[mid]) {
     90:         break;
     91:       } else if (idex < petsc[mid]) {
     92:         high = mid - 1;
     93:       } else {
     94:         low  = mid + 1;
     95:       }
     96:     }
     97:     if (low > high) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Invalid input index %D", idex);
     98:     ia[i] = app[perm[mid]];
     99:   }
    100:   return(0);
    101: }

    105: PetscErrorCode AOApplicationToPetsc_Mapping(AO ao, PetscInt n, PetscInt *ia)
    106: {
    107:   AO_Mapping *aomap = (AO_Mapping *) ao->data;
    108:   PetscInt   *app   = aomap->app;
    109:   PetscInt   *petsc = aomap->petsc;
    110:   PetscInt   *perm  = aomap->appPerm;
    111:   PetscInt   N     = aomap->N;
    112:   PetscInt   low, high, mid=0;
    113:   PetscInt   idex;
    114:   PetscInt   i;

    116:   /* It would be possible to use a single bisection search, which
    117:      recursively divided the indices to be converted, and searched
    118:      partitions which contained an index. This would result in
    119:      better running times if indices are clustered.
    120:   */
    122:   for(i = 0; i < n; i++) {
    123:     idex = ia[i];
    124:     if (idex < 0) continue;
    125:     /* Use bisection since the array is sorted */
    126:     low  = 0;
    127:     high = N - 1;
    128:     while (low <= high) {
    129:       mid = (low + high)/2;
    130:       if (idex == app[mid]) {
    131:         break;
    132:       } else if (idex < app[mid]) {
    133:         high = mid - 1;
    134:       } else {
    135:         low  = mid + 1;
    136:       }
    137:     }
    138:     if (low > high) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Invalid input index %D", idex);
    139:     ia[i] = petsc[perm[mid]];
    140:   }
    141:   return(0);
    142: }

    144: static struct _AOOps AOps = {AOView_Mapping,
    145:                              AODestroy_Mapping,
    146:                              AOPetscToApplication_Mapping,
    147:                              AOApplicationToPetsc_Mapping,
    148:                              PETSC_NULL,
    149:                              PETSC_NULL,
    150:                              PETSC_NULL,
    151:                              PETSC_NULL};

    155: /*@C
    156:   AOMappingHasApplicationIndex - Searches for the supplied application index.

    158:   Input Parameters:
    159: + ao       - The AOMapping
    160: - index    - The application index

    162:   Output Parameter:
    163: . hasIndex - Flag is PETSC_TRUE if the index exists

    165:   Level: intermediate

    167: .keywords: AO, index
    168: .seealso: AOMappingHasPetscIndex(), AOCreateMapping()
    169: @*/
    170: PetscErrorCode PETSCDM_DLLEXPORT AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscTruth *hasIndex)
    171: {
    172:   AO_Mapping *aomap;
    173:   PetscInt   *app;
    174:   PetscInt   low, high, mid;

    179:   aomap = (AO_Mapping *) ao->data;
    180:   app   = aomap->app;
    181:   /* Use bisection since the array is sorted */
    182:   low  = 0;
    183:   high = aomap->N - 1;
    184:   while (low <= high) {
    185:     mid = (low + high)/2;
    186:     if (idex == app[mid]) {
    187:       break;
    188:     } else if (idex < app[mid]) {
    189:       high = mid - 1;
    190:     } else {
    191:       low  = mid + 1;
    192:     }
    193:   }
    194:   if (low > high) {
    195:     *hasIndex = PETSC_FALSE;
    196:   } else {
    197:     *hasIndex = PETSC_TRUE;
    198:   }
    199:   return(0);
    200: }

    204: /*@C
    205:   AOMappingHasPetscIndex - Searches for the supplied petsc index.

    207:   Input Parameters:
    208: + ao       - The AOMapping
    209: - index    - The petsc index

    211:   Output Parameter:
    212: . hasIndex - Flag is PETSC_TRUE if the index exists

    214:   Level: intermediate

    216: .keywords: AO, index
    217: .seealso: AOMappingHasApplicationIndex(), AOCreateMapping()
    218: @*/
    219: PetscErrorCode PETSCDM_DLLEXPORT AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscTruth *hasIndex)
    220: {
    221:   AO_Mapping *aomap;
    222:   PetscInt   *petsc;
    223:   PetscInt   low, high, mid;

    228:   aomap = (AO_Mapping *) ao->data;
    229:   petsc = aomap->petsc;
    230:   /* Use bisection since the array is sorted */
    231:   low  = 0;
    232:   high = aomap->N - 1;
    233:   while (low <= high) {
    234:     mid = (low + high)/2;
    235:     if (idex == petsc[mid]) {
    236:       break;
    237:     } else if (idex < petsc[mid]) {
    238:       high = mid - 1;
    239:     } else {
    240:       low  = mid + 1;
    241:     }
    242:   }
    243:   if (low > high) {
    244:     *hasIndex = PETSC_FALSE;
    245:   } else {
    246:     *hasIndex = PETSC_TRUE;
    247:   }
    248:   return(0);
    249: }

    253: /*@C
    254:   AOCreateMapping - Creates a basic application mapping using two integer arrays.

    256:   Input Parameters:
    257: + comm    - MPI communicator that is to share AO
    258: . napp    - size of integer arrays
    259: . myapp   - integer array that defines an ordering
    260: - mypetsc - integer array that defines another ordering

    262:   Output Parameter:
    263: . aoout   - the new application mapping

    265:   Options Database Key:
    266: $ -ao_view : call AOView() at the conclusion of AOCreateMapping()

    268:   Level: beginner

    270: .keywords: AO, create
    271: .seealso: AOCreateDebug(), AOCreateBasic(), AOCreateMappingIS(), AODestroy()
    272: @*/
    273: PetscErrorCode PETSCDM_DLLEXPORT AOCreateMapping(MPI_Comm comm,PetscInt napp,const PetscInt myapp[],const PetscInt mypetsc[],AO *aoout)
    274: {
    275:   AO             ao;
    276:   AO_Mapping     *aomap;
    277:   PetscInt       *allpetsc,  *allapp;
    278:   PetscInt       *petscPerm, *appPerm;
    279:   PetscInt       *petsc;
    280:   PetscMPIInt    size, rank,*lens, *disp,nnapp;
    281:   PetscInt       N, start;
    282:   PetscInt       i;
    283:   PetscTruth     opt;

    288:   *aoout = 0;
    289: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
    290:   DMInitializePackage(PETSC_NULL);
    291: #endif

    293:   PetscHeaderCreate(ao, _p_AO, struct _AOOps, AO_COOKIE, AO_MAPPING, "AO", comm, AODestroy, AOView);
    294:   PetscNew(AO_Mapping, &aomap);
    295:   PetscLogObjectMemory(ao, sizeof(struct _p_AO) + sizeof(AO_Mapping));
    296:   PetscMemcpy(ao->ops, &AOps, sizeof(AOps));
    297:   ao->data = (void*) aomap;

    299:   /* transmit all lengths to all processors */
    300:   MPI_Comm_size(comm, &size);
    301:   MPI_Comm_rank(comm, &rank);
    302:   PetscMalloc(2*size * sizeof(PetscMPIInt), &lens);
    303:   disp  = lens + size;
    304:   nnapp = napp;
    305:   MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm);
    306:   N    = 0;
    307:   for(i = 0; i < size; i++) {
    308:     disp[i] = N;
    309:     N += lens[i];
    310:   }
    311:   aomap->N = N;
    312:   ao->N    = N;
    313:   ao->n    = N;

    315:   /* If mypetsc is 0 then use "natural" numbering */
    316:   if (!mypetsc) {
    317:     start = disp[rank];
    318:     PetscMalloc((napp+1) * sizeof(PetscInt), &petsc);
    319:     for(i = 0; i < napp; i++) {
    320:       petsc[i] = start + i;
    321:     }
    322:   } else {
    323:     petsc = (PetscInt*)mypetsc;
    324:   }

    326:   /* get all indices on all processors */
    327:   PetscMalloc(N*4 * sizeof(PetscInt), &allapp);
    328:   appPerm   = allapp   + N;
    329:   allpetsc  = appPerm  + N;
    330:   petscPerm = allpetsc + N;
    331:   MPI_Allgatherv((void*)myapp, napp, MPIU_INT, allapp,   lens, disp, MPIU_INT, comm);
    332:   MPI_Allgatherv((void*)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm);
    333:   PetscFree(lens);

    335:   /* generate a list of application and PETSc node numbers */
    336:   PetscMalloc(N*4 * sizeof(PetscInt), &aomap->app);
    337:   PetscLogObjectMemory(ao, 4*N * sizeof(PetscInt));
    338:   aomap->appPerm   = aomap->app     + N;
    339:   aomap->petsc     = aomap->appPerm + N;
    340:   aomap->petscPerm = aomap->petsc   + N;
    341:   for(i = 0; i < N; i++) {
    342:     appPerm[i]   = i;
    343:     petscPerm[i] = i;
    344:   }
    345:   PetscSortIntWithPermutation(N, allpetsc, petscPerm);
    346:   PetscSortIntWithPermutation(N, allapp,   appPerm);
    347:   /* Form sorted arrays of indices */
    348:   for(i = 0; i < N; i++) {
    349:     aomap->app[i]   = allapp[appPerm[i]];
    350:     aomap->petsc[i] = allpetsc[petscPerm[i]];
    351:   }
    352:   /* Invert petscPerm[] into aomap->petscPerm[] */
    353:   for(i = 0; i < N; i++) {
    354:     aomap->petscPerm[petscPerm[i]] = i;
    355:   }
    356:   /* Form map between aomap->app[] and aomap->petsc[] */
    357:   for(i = 0; i < N; i++) {
    358:     aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];
    359:   }
    360:   /* Invert appPerm[] into allapp[] */
    361:   for(i = 0; i < N; i++) {
    362:     allapp[appPerm[i]] = i;
    363:   }
    364:   /* Form map between aomap->petsc[] and aomap->app[] */
    365:   for(i = 0; i < N; i++) {
    366:     aomap->petscPerm[i] = allapp[petscPerm[i]];
    367:   }
    368: #ifdef PETSC_USE_DEBUG
    369:   /* Check that the permutations are complementary */
    370:   for(i = 0; i < N; i++) {
    371:     if (i != aomap->appPerm[aomap->petscPerm[i]])
    372:       SETERRQ(PETSC_ERR_PLIB, "Invalid ordering");
    373:   }
    374: #endif
    375:   /* Cleanup */
    376:   if (!mypetsc) {
    377:     PetscFree(petsc);
    378:   }
    379:   PetscFree(allapp);

    381:   PetscOptionsHasName(PETSC_NULL, "-ao_view", &opt);
    382:   if (opt) {
    383:     AOView(ao, PETSC_VIEWER_STDOUT_SELF);
    384:   }

    386:   *aoout = ao;
    387:   return(0);
    388: }

    392: /*@C
    393:   AOCreateMappingIS - Creates a basic application ordering using two index sets.

    395:   Input Parameters:
    396: + comm    - MPI communicator that is to share AO
    397: . isapp   - index set that defines an ordering
    398: - ispetsc - index set that defines another ordering

    400:   Output Parameter:
    401: . aoout   - the new application ordering

    403:   Options Database Key:
    404: $ -ao_view : call AOView() at the conclusion of AOCreateMappingIS()

    406:   Level: beginner

    408: .keywords: AO, create
    409: .seealso: AOCreateBasic(), AOCreateMapping(), AODestroy()
    410: @*/
    411: PetscErrorCode PETSCDM_DLLEXPORT AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
    412: {
    413:   MPI_Comm       comm;
    414:   PetscInt       *mypetsc, *myapp;
    415:   PetscInt       napp, npetsc;

    419:   PetscObjectGetComm((PetscObject) isapp, &comm);
    420:   ISGetSize(isapp, &napp);
    421:   if (ispetsc) {
    422:     ISGetSize(ispetsc, &npetsc);
    423:     if (napp != npetsc) SETERRQ(PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
    424:     ISGetIndices(ispetsc, &mypetsc);
    425:   } else {
    426:     mypetsc = NULL;
    427:   }
    428:   ISGetIndices(isapp, &myapp);

    430:   AOCreateMapping(comm, napp, myapp, mypetsc, aoout);

    432:   ISRestoreIndices(isapp, &myapp);
    433:   if (ispetsc) {
    434:     ISRestoreIndices(ispetsc, &mypetsc);
    435:   }
    436:   return(0);
    437: }