Actual source code: aomapping.c

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

      7: #include <../src/vec/is/ao/aoimpl.h>

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

     17: static PetscErrorCode AODestroy_Mapping(AO ao)
     18: {
     19:   AO_Mapping *aomap = (AO_Mapping *)ao->data;

     21:   PetscFunctionBegin;
     22:   PetscCall(PetscFree4(aomap->app, aomap->appPerm, aomap->petsc, aomap->petscPerm));
     23:   PetscCall(PetscFree(aomap));
     24:   PetscFunctionReturn(PETSC_SUCCESS);
     25: }

     27: static PetscErrorCode AOView_Mapping(AO ao, PetscViewer viewer)
     28: {
     29:   AO_Mapping *aomap = (AO_Mapping *)ao->data;
     30:   PetscMPIInt rank;
     31:   PetscInt    i;
     32:   PetscBool   iascii;

     34:   PetscFunctionBegin;
     35:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank));
     36:   if (rank) PetscFunctionReturn(PETSC_SUCCESS);
     37:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
     38:   if (iascii) {
     39:     PetscCall(PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %" PetscInt_FMT "\n", aomap->N));
     40:     PetscCall(PetscViewerASCIIPrintf(viewer, "   App.   PETSc\n"));
     41:     for (i = 0; i < aomap->N; i++) PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT "   %" PetscInt_FMT "    %" PetscInt_FMT "\n", i, aomap->app[i], aomap->petsc[aomap->appPerm[i]]));
     42:   }
     43:   PetscFunctionReturn(PETSC_SUCCESS);
     44: }

     46: static PetscErrorCode AOPetscToApplication_Mapping(AO ao, PetscInt n, PetscInt *ia)
     47: {
     48:   AO_Mapping *aomap = (AO_Mapping *)ao->data;
     49:   PetscInt   *app   = aomap->app;
     50:   PetscInt   *petsc = aomap->petsc;
     51:   PetscInt   *perm  = aomap->petscPerm;
     52:   PetscInt    N     = aomap->N;
     53:   PetscInt    low, high, mid = 0;
     54:   PetscInt    idex;
     55:   PetscInt    i;

     57:   /* It would be possible to use a single bisection search, which
     58:      recursively divided the indices to be converted, and searched
     59:      partitions which contained an index. This would result in
     60:      better running times if indices are clustered.
     61:   */
     62:   PetscFunctionBegin;
     63:   for (i = 0; i < n; i++) {
     64:     idex = ia[i];
     65:     if (idex < 0) continue;
     66:     /* Use bisection since the array is sorted */
     67:     low  = 0;
     68:     high = N - 1;
     69:     while (low <= high) {
     70:       mid = (low + high) / 2;
     71:       if (idex == petsc[mid]) break;
     72:       else if (idex < petsc[mid]) high = mid - 1;
     73:       else low = mid + 1;
     74:     }
     75:     if (low > high) ia[i] = -1;
     76:     else ia[i] = app[perm[mid]];
     77:   }
     78:   PetscFunctionReturn(PETSC_SUCCESS);
     79: }

     81: static PetscErrorCode AOApplicationToPetsc_Mapping(AO ao, PetscInt n, PetscInt *ia)
     82: {
     83:   AO_Mapping *aomap = (AO_Mapping *)ao->data;
     84:   PetscInt   *app   = aomap->app;
     85:   PetscInt   *petsc = aomap->petsc;
     86:   PetscInt   *perm  = aomap->appPerm;
     87:   PetscInt    N     = aomap->N;
     88:   PetscInt    low, high, mid = 0;
     89:   PetscInt    idex;
     90:   PetscInt    i;

     92:   /* It would be possible to use a single bisection search, which
     93:      recursively divided the indices to be converted, and searched
     94:      partitions which contained an index. This would result in
     95:      better running times if indices are clustered.
     96:   */
     97:   PetscFunctionBegin;
     98:   for (i = 0; i < n; i++) {
     99:     idex = ia[i];
    100:     if (idex < 0) continue;
    101:     /* Use bisection since the array is sorted */
    102:     low  = 0;
    103:     high = N - 1;
    104:     while (low <= high) {
    105:       mid = (low + high) / 2;
    106:       if (idex == app[mid]) break;
    107:       else if (idex < app[mid]) high = mid - 1;
    108:       else low = mid + 1;
    109:     }
    110:     if (low > high) ia[i] = -1;
    111:     else ia[i] = petsc[perm[mid]];
    112:   }
    113:   PetscFunctionReturn(PETSC_SUCCESS);
    114: }

    116: static const struct _AOOps AOps = {
    117:   PetscDesignatedInitializer(view, AOView_Mapping),
    118:   PetscDesignatedInitializer(destroy, AODestroy_Mapping),
    119:   PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_Mapping),
    120:   PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_Mapping),
    121:   PetscDesignatedInitializer(petsctoapplicationpermuteint, NULL),
    122:   PetscDesignatedInitializer(applicationtopetscpermuteint, NULL),
    123:   PetscDesignatedInitializer(petsctoapplicationpermutereal, NULL),
    124:   PetscDesignatedInitializer(applicationtopetscpermutereal, NULL),
    125: };

    127: /*@
    128:   AOMappingHasApplicationIndex - Checks if an `AO` has a requested application index.

    130:   Not Collective

    132:   Input Parameters:
    133: + ao   - The `AO`
    134: - idex - The application index

    136:   Output Parameter:
    137: . hasIndex - Flag is `PETSC_TRUE` if the index exists

    139:   Level: intermediate

    141:   Developer Notes:
    142:   The name of the function is wrong, it should be `AOHasApplicationIndex`

    144: .seealso: [](sec_ao), `AOMappingHasPetscIndex()`, `AOCreateMapping()`, `AO`
    145: @*/
    146: PetscErrorCode AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
    147: {
    148:   AO_Mapping *aomap;
    149:   PetscInt   *app;
    150:   PetscInt    low, high, mid;

    152:   PetscFunctionBegin;
    154:   PetscAssertPointer(hasIndex, 3);
    155:   aomap = (AO_Mapping *)ao->data;
    156:   app   = aomap->app;
    157:   /* Use bisection since the array is sorted */
    158:   low  = 0;
    159:   high = aomap->N - 1;
    160:   while (low <= high) {
    161:     mid = (low + high) / 2;
    162:     if (idex == app[mid]) break;
    163:     else if (idex < app[mid]) high = mid - 1;
    164:     else low = mid + 1;
    165:   }
    166:   if (low > high) *hasIndex = PETSC_FALSE;
    167:   else *hasIndex = PETSC_TRUE;
    168:   PetscFunctionReturn(PETSC_SUCCESS);
    169: }

    171: /*@
    172:   AOMappingHasPetscIndex - checks if an `AO` has a requested petsc index.

    174:   Not Collective

    176:   Input Parameters:
    177: + ao   - The `AO`
    178: - idex - The petsc index

    180:   Output Parameter:
    181: . hasIndex - Flag is `PETSC_TRUE` if the index exists

    183:   Level: intermediate

    185:   Developer Notes:
    186:   The name of the function is wrong, it should be `AOHasPetscIndex`

    188: .seealso: [](sec_ao), `AOMappingHasApplicationIndex()`, `AOCreateMapping()`
    189: @*/
    190: PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool *hasIndex)
    191: {
    192:   AO_Mapping *aomap;
    193:   PetscInt   *petsc;
    194:   PetscInt    low, high, mid;

    196:   PetscFunctionBegin;
    198:   PetscAssertPointer(hasIndex, 3);
    199:   aomap = (AO_Mapping *)ao->data;
    200:   petsc = aomap->petsc;
    201:   /* Use bisection since the array is sorted */
    202:   low  = 0;
    203:   high = aomap->N - 1;
    204:   while (low <= high) {
    205:     mid = (low + high) / 2;
    206:     if (idex == petsc[mid]) break;
    207:     else if (idex < petsc[mid]) high = mid - 1;
    208:     else low = mid + 1;
    209:   }
    210:   if (low > high) *hasIndex = PETSC_FALSE;
    211:   else *hasIndex = PETSC_TRUE;
    212:   PetscFunctionReturn(PETSC_SUCCESS);
    213: }

    215: /*@
    216:   AOCreateMapping - Creates an application mapping using two integer arrays.

    218:   Input Parameters:
    219: + comm    - MPI communicator that is to share the `AO`
    220: . napp    - size of integer arrays
    221: . myapp   - integer array that defines an ordering
    222: - mypetsc - integer array that defines another ordering (may be `NULL` to indicate the identity ordering)

    224:   Output Parameter:
    225: . aoout - the new application mapping

    227:   Options Database Key:
    228: . -ao_view - call `AOView()` at the conclusion of `AOCreateMapping()`

    230:   Level: beginner

    232:   Note:
    233:   The arrays `myapp` and `mypetsc` need NOT contain the all the integers 0 to `napp`-1, that is there CAN be "holes"  in the indices.
    234:   Use `AOCreateBasic()` or `AOCreateBasicIS()` if they do not have holes for better performance.

    236: .seealso: [](sec_ao), `AOCreateBasic()`, `AOCreateMappingIS()`, `AODestroy()`
    237: @*/
    238: PetscErrorCode AOCreateMapping(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout)
    239: {
    240:   AO          ao;
    241:   AO_Mapping *aomap;
    242:   PetscInt   *allpetsc, *allapp;
    243:   PetscInt   *petscPerm, *appPerm;
    244:   PetscInt   *petsc;
    245:   PetscMPIInt size, rank, *lens, *disp, nnapp;
    246:   PetscCount  N;
    247:   PetscInt    start;

    249:   PetscFunctionBegin;
    250:   PetscAssertPointer(aoout, 5);
    251:   *aoout = NULL;
    252:   PetscCall(AOInitializePackage());

    254:   PetscCall(PetscHeaderCreate(ao, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView));
    255:   PetscCall(PetscNew(&aomap));
    256:   ao->ops[0] = AOps;
    257:   ao->data   = (void *)aomap;

    259:   /* transmit all lengths to all processors */
    260:   PetscCallMPI(MPI_Comm_size(comm, &size));
    261:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    262:   PetscCall(PetscMalloc2(size, &lens, size, &disp));
    263:   PetscCall(PetscMPIIntCast(napp, &nnapp));
    264:   PetscCallMPI(MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm));
    265:   N = 0;
    266:   for (PetscMPIInt i = 0; i < size; i++) {
    267:     PetscCall(PetscMPIIntCast(N, &disp[i]));
    268:     N += lens[i];
    269:   }
    270:   PetscCall(PetscIntCast(N, &aomap->N));
    271:   ao->N = (PetscInt)N;
    272:   ao->n = (PetscInt)N;

    274:   /* If mypetsc is 0 then use "natural" numbering */
    275:   if (!mypetsc) {
    276:     start = disp[rank];
    277:     PetscCall(PetscMalloc1(napp + 1, &petsc));
    278:     for (PetscInt i = 0; i < napp; i++) petsc[i] = start + i;
    279:   } else {
    280:     petsc = (PetscInt *)mypetsc;
    281:   }

    283:   /* get all indices on all processors */
    284:   PetscCall(PetscMalloc4(N, &allapp, N, &appPerm, N, &allpetsc, N, &petscPerm));
    285:   PetscCallMPI(MPI_Allgatherv((void *)myapp, nnapp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm));
    286:   PetscCallMPI(MPI_Allgatherv((void *)petsc, nnapp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm));
    287:   PetscCall(PetscFree2(lens, disp));

    289:   /* generate a list of application and PETSc node numbers */
    290:   PetscCall(PetscMalloc4(N, &aomap->app, N, &aomap->appPerm, N, &aomap->petsc, N, &aomap->petscPerm));
    291:   for (PetscInt i = 0; i < N; i++) {
    292:     appPerm[i]   = i;
    293:     petscPerm[i] = i;
    294:   }
    295:   PetscCall(PetscSortIntWithPermutation((PetscInt)N, allpetsc, petscPerm));
    296:   PetscCall(PetscSortIntWithPermutation((PetscInt)N, allapp, appPerm));
    297:   /* Form sorted arrays of indices */
    298:   for (PetscInt i = 0; i < N; i++) {
    299:     aomap->app[i]   = allapp[appPerm[i]];
    300:     aomap->petsc[i] = allpetsc[petscPerm[i]];
    301:   }
    302:   /* Invert petscPerm[] into aomap->petscPerm[] */
    303:   for (PetscInt i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i;

    305:   /* Form map between aomap->app[] and aomap->petsc[] */
    306:   for (PetscInt i = 0; i < N; i++) aomap->appPerm[i] = aomap->petscPerm[appPerm[i]];

    308:   /* Invert appPerm[] into allapp[] */
    309:   for (PetscInt i = 0; i < N; i++) allapp[appPerm[i]] = i;

    311:   /* Form map between aomap->petsc[] and aomap->app[] */
    312:   for (PetscInt i = 0; i < N; i++) aomap->petscPerm[i] = allapp[petscPerm[i]];

    314:   if (PetscDefined(USE_DEBUG)) {
    315:     /* Check that the permutations are complementary */
    316:     for (PetscInt i = 0; i < N; i++) PetscCheck(i == aomap->appPerm[aomap->petscPerm[i]], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid ordering");
    317:   }
    318:   /* Cleanup */
    319:   if (!mypetsc) PetscCall(PetscFree(petsc));
    320:   PetscCall(PetscFree4(allapp, appPerm, allpetsc, petscPerm));
    321:   PetscCall(AOViewFromOptions(ao, NULL, "-ao_view"));
    322:   *aoout = ao;
    323:   PetscFunctionReturn(PETSC_SUCCESS);
    324: }

    326: /*@
    327:   AOCreateMappingIS - Creates an application mapping using two index sets.

    329:   Input Parameters:
    330: + isapp   - index set that defines an ordering
    331: - ispetsc - index set that defines another ordering, maybe `NULL` for identity `IS`

    333:   Output Parameter:
    334: . aoout - the new application ordering

    336:   Options Database Key:
    337: . -ao_view - call `AOView()` at the conclusion of `AOCreateMappingIS()`

    339:   Level: beginner

    341:   Note:
    342:   The index sets `isapp` and `ispetsc` need NOT contain the all the integers 0 to N-1, that is there CAN be "holes"  in the indices.
    343:   Use `AOCreateBasic()` or `AOCreateBasicIS()` if they do not have holes for better performance.

    345: .seealso: [](sec_ao), [](sec_scatter), `AOCreateBasic()`, `AOCreateMapping()`, `AODestroy()`
    346: @*/
    347: PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout)
    348: {
    349:   MPI_Comm        comm;
    350:   const PetscInt *mypetsc, *myapp;
    351:   PetscInt        napp, npetsc;

    353:   PetscFunctionBegin;
    354:   PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm));
    355:   PetscCall(ISGetLocalSize(isapp, &napp));
    356:   if (ispetsc) {
    357:     PetscCall(ISGetLocalSize(ispetsc, &npetsc));
    358:     PetscCheck(napp == npetsc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local IS lengths must match");
    359:     PetscCall(ISGetIndices(ispetsc, &mypetsc));
    360:   } else {
    361:     mypetsc = NULL;
    362:   }
    363:   PetscCall(ISGetIndices(isapp, &myapp));

    365:   PetscCall(AOCreateMapping(comm, napp, myapp, mypetsc, aoout));

    367:   PetscCall(ISRestoreIndices(isapp, &myapp));
    368:   if (ispetsc) PetscCall(ISRestoreIndices(ispetsc, &mypetsc));
    369:   PetscFunctionReturn(PETSC_SUCCESS);
    370: }