1 2 /* 3 These AO application ordering routines do not require that the input 4 be a permutation, but merely a 1-1 mapping. This implementation still 5 keeps the entire ordering on each processor. 6 */ 7 8 #include <../src/vec/is/ao/aoimpl.h> /*I "petscao.h" I*/ 9 10 typedef struct { 11 PetscInt N; 12 PetscInt *app; /* app[i] is the partner for petsc[appPerm[i]] */ 13 PetscInt *appPerm; 14 PetscInt *petsc; /* petsc[j] is the partner for app[petscPerm[j]] */ 15 PetscInt *petscPerm; 16 } AO_Mapping; 17 18 PetscErrorCode AODestroy_Mapping(AO ao) 19 { 20 AO_Mapping *aomap = (AO_Mapping *)ao->data; 21 22 PetscFunctionBegin; 23 PetscCall(PetscFree4(aomap->app, aomap->appPerm, aomap->petsc, aomap->petscPerm)); 24 PetscCall(PetscFree(aomap)); 25 PetscFunctionReturn(PETSC_SUCCESS); 26 } 27 28 PetscErrorCode AOView_Mapping(AO ao, PetscViewer viewer) 29 { 30 AO_Mapping *aomap = (AO_Mapping *)ao->data; 31 PetscMPIInt rank; 32 PetscInt i; 33 PetscBool iascii; 34 35 PetscFunctionBegin; 36 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank)); 37 if (rank) PetscFunctionReturn(PETSC_SUCCESS); 38 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 39 if (iascii) { 40 PetscCall(PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %" PetscInt_FMT "\n", aomap->N)); 41 PetscCall(PetscViewerASCIIPrintf(viewer, " App. PETSc\n")); 42 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]])); 43 } 44 PetscFunctionReturn(PETSC_SUCCESS); 45 } 46 47 PetscErrorCode AOPetscToApplication_Mapping(AO ao, PetscInt n, PetscInt *ia) 48 { 49 AO_Mapping *aomap = (AO_Mapping *)ao->data; 50 PetscInt *app = aomap->app; 51 PetscInt *petsc = aomap->petsc; 52 PetscInt *perm = aomap->petscPerm; 53 PetscInt N = aomap->N; 54 PetscInt low, high, mid = 0; 55 PetscInt idex; 56 PetscInt i; 57 58 /* It would be possible to use a single bisection search, which 59 recursively divided the indices to be converted, and searched 60 partitions which contained an index. This would result in 61 better running times if indices are clustered. 62 */ 63 PetscFunctionBegin; 64 for (i = 0; i < n; i++) { 65 idex = ia[i]; 66 if (idex < 0) continue; 67 /* Use bisection since the array is sorted */ 68 low = 0; 69 high = N - 1; 70 while (low <= high) { 71 mid = (low + high) / 2; 72 if (idex == petsc[mid]) break; 73 else if (idex < petsc[mid]) high = mid - 1; 74 else low = mid + 1; 75 } 76 if (low > high) ia[i] = -1; 77 else ia[i] = app[perm[mid]]; 78 } 79 PetscFunctionReturn(PETSC_SUCCESS); 80 } 81 82 PetscErrorCode AOApplicationToPetsc_Mapping(AO ao, PetscInt n, PetscInt *ia) 83 { 84 AO_Mapping *aomap = (AO_Mapping *)ao->data; 85 PetscInt *app = aomap->app; 86 PetscInt *petsc = aomap->petsc; 87 PetscInt *perm = aomap->appPerm; 88 PetscInt N = aomap->N; 89 PetscInt low, high, mid = 0; 90 PetscInt idex; 91 PetscInt i; 92 93 /* It would be possible to use a single bisection search, which 94 recursively divided the indices to be converted, and searched 95 partitions which contained an index. This would result in 96 better running times if indices are clustered. 97 */ 98 PetscFunctionBegin; 99 for (i = 0; i < n; i++) { 100 idex = ia[i]; 101 if (idex < 0) continue; 102 /* Use bisection since the array is sorted */ 103 low = 0; 104 high = N - 1; 105 while (low <= high) { 106 mid = (low + high) / 2; 107 if (idex == app[mid]) break; 108 else if (idex < app[mid]) high = mid - 1; 109 else low = mid + 1; 110 } 111 if (low > high) ia[i] = -1; 112 else ia[i] = petsc[perm[mid]]; 113 } 114 PetscFunctionReturn(PETSC_SUCCESS); 115 } 116 117 static struct _AOOps AOps = { 118 PetscDesignatedInitializer(view, AOView_Mapping), 119 PetscDesignatedInitializer(destroy, AODestroy_Mapping), 120 PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_Mapping), 121 PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_Mapping), 122 }; 123 124 /*@C 125 AOMappingHasApplicationIndex - Checks if an `AO` has a requested application index. 126 127 Not Collective 128 129 Input Parameters: 130 + ao - The `AO` 131 - idex - The application index 132 133 Output Parameter: 134 . hasIndex - Flag is `PETSC_TRUE` if the index exists 135 136 Level: intermediate 137 138 Developer Notes: 139 The name of the function is wrong, it should be `AOHasApplicationIndex` 140 141 .seealso: [](sec_ao), `AOMappingHasPetscIndex()`, `AOCreateMapping()`, `AO` 142 @*/ 143 PetscErrorCode AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool *hasIndex) 144 { 145 AO_Mapping *aomap; 146 PetscInt *app; 147 PetscInt low, high, mid; 148 149 PetscFunctionBegin; 150 PetscValidHeaderSpecific(ao, AO_CLASSID, 1); 151 PetscAssertPointer(hasIndex, 3); 152 aomap = (AO_Mapping *)ao->data; 153 app = aomap->app; 154 /* Use bisection since the array is sorted */ 155 low = 0; 156 high = aomap->N - 1; 157 while (low <= high) { 158 mid = (low + high) / 2; 159 if (idex == app[mid]) break; 160 else if (idex < app[mid]) high = mid - 1; 161 else low = mid + 1; 162 } 163 if (low > high) *hasIndex = PETSC_FALSE; 164 else *hasIndex = PETSC_TRUE; 165 PetscFunctionReturn(PETSC_SUCCESS); 166 } 167 168 /*@ 169 AOMappingHasPetscIndex - checks if an `AO` has a requested petsc index. 170 171 Not Collective 172 173 Input Parameters: 174 + ao - The `AO` 175 - idex - The petsc index 176 177 Output Parameter: 178 . hasIndex - Flag is `PETSC_TRUE` if the index exists 179 180 Level: intermediate 181 182 Developer Notes: 183 The name of the function is wrong, it should be `AOHasPetscIndex` 184 185 .seealso: [](sec_ao), `AOMappingHasApplicationIndex()`, `AOCreateMapping()` 186 @*/ 187 PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool *hasIndex) 188 { 189 AO_Mapping *aomap; 190 PetscInt *petsc; 191 PetscInt low, high, mid; 192 193 PetscFunctionBegin; 194 PetscValidHeaderSpecific(ao, AO_CLASSID, 1); 195 PetscAssertPointer(hasIndex, 3); 196 aomap = (AO_Mapping *)ao->data; 197 petsc = aomap->petsc; 198 /* Use bisection since the array is sorted */ 199 low = 0; 200 high = aomap->N - 1; 201 while (low <= high) { 202 mid = (low + high) / 2; 203 if (idex == petsc[mid]) break; 204 else if (idex < petsc[mid]) high = mid - 1; 205 else low = mid + 1; 206 } 207 if (low > high) *hasIndex = PETSC_FALSE; 208 else *hasIndex = PETSC_TRUE; 209 PetscFunctionReturn(PETSC_SUCCESS); 210 } 211 212 /*@C 213 AOCreateMapping - Creates an application mapping using two integer arrays. 214 215 Input Parameters: 216 + comm - MPI communicator that is to share the `AO` 217 . napp - size of integer arrays 218 . myapp - integer array that defines an ordering 219 - mypetsc - integer array that defines another ordering (may be NULL to indicate the identity ordering) 220 221 Output Parameter: 222 . aoout - the new application mapping 223 224 Options Database Key: 225 . -ao_view - call `AOView()` at the conclusion of `AOCreateMapping()` 226 227 Level: beginner 228 229 Note: 230 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. 231 Use `AOCreateBasic()` or `AOCreateBasicIS()` if they do not have holes for better performance. 232 233 .seealso: [](sec_ao), `AOCreateBasic()`, `AOCreateMappingIS()`, `AODestroy()` 234 @*/ 235 PetscErrorCode AOCreateMapping(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout) 236 { 237 AO ao; 238 AO_Mapping *aomap; 239 PetscInt *allpetsc, *allapp; 240 PetscInt *petscPerm, *appPerm; 241 PetscInt *petsc; 242 PetscMPIInt size, rank, *lens, *disp, nnapp; 243 PetscInt N, start; 244 PetscInt i; 245 246 PetscFunctionBegin; 247 PetscAssertPointer(aoout, 5); 248 *aoout = NULL; 249 PetscCall(AOInitializePackage()); 250 251 PetscCall(PetscHeaderCreate(ao, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView)); 252 PetscCall(PetscNew(&aomap)); 253 ao->ops[0] = AOps; 254 ao->data = (void *)aomap; 255 256 /* transmit all lengths to all processors */ 257 PetscCallMPI(MPI_Comm_size(comm, &size)); 258 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 259 PetscCall(PetscMalloc2(size, &lens, size, &disp)); 260 nnapp = napp; 261 PetscCallMPI(MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm)); 262 N = 0; 263 for (i = 0; i < size; i++) { 264 disp[i] = N; 265 N += lens[i]; 266 } 267 aomap->N = N; 268 ao->N = N; 269 ao->n = N; 270 271 /* If mypetsc is 0 then use "natural" numbering */ 272 if (!mypetsc) { 273 start = disp[rank]; 274 PetscCall(PetscMalloc1(napp + 1, &petsc)); 275 for (i = 0; i < napp; i++) petsc[i] = start + i; 276 } else { 277 petsc = (PetscInt *)mypetsc; 278 } 279 280 /* get all indices on all processors */ 281 PetscCall(PetscMalloc4(N, &allapp, N, &appPerm, N, &allpetsc, N, &petscPerm)); 282 PetscCallMPI(MPI_Allgatherv((void *)myapp, napp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm)); 283 PetscCallMPI(MPI_Allgatherv((void *)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm)); 284 PetscCall(PetscFree2(lens, disp)); 285 286 /* generate a list of application and PETSc node numbers */ 287 PetscCall(PetscMalloc4(N, &aomap->app, N, &aomap->appPerm, N, &aomap->petsc, N, &aomap->petscPerm)); 288 for (i = 0; i < N; i++) { 289 appPerm[i] = i; 290 petscPerm[i] = i; 291 } 292 PetscCall(PetscSortIntWithPermutation(N, allpetsc, petscPerm)); 293 PetscCall(PetscSortIntWithPermutation(N, allapp, appPerm)); 294 /* Form sorted arrays of indices */ 295 for (i = 0; i < N; i++) { 296 aomap->app[i] = allapp[appPerm[i]]; 297 aomap->petsc[i] = allpetsc[petscPerm[i]]; 298 } 299 /* Invert petscPerm[] into aomap->petscPerm[] */ 300 for (i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i; 301 302 /* Form map between aomap->app[] and aomap->petsc[] */ 303 for (i = 0; i < N; i++) aomap->appPerm[i] = aomap->petscPerm[appPerm[i]]; 304 305 /* Invert appPerm[] into allapp[] */ 306 for (i = 0; i < N; i++) allapp[appPerm[i]] = i; 307 308 /* Form map between aomap->petsc[] and aomap->app[] */ 309 for (i = 0; i < N; i++) aomap->petscPerm[i] = allapp[petscPerm[i]]; 310 311 if (PetscDefined(USE_DEBUG)) { 312 /* Check that the permutations are complementary */ 313 for (i = 0; i < N; i++) PetscCheck(i == aomap->appPerm[aomap->petscPerm[i]], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid ordering"); 314 } 315 /* Cleanup */ 316 if (!mypetsc) PetscCall(PetscFree(petsc)); 317 PetscCall(PetscFree4(allapp, appPerm, allpetsc, petscPerm)); 318 319 PetscCall(AOViewFromOptions(ao, NULL, "-ao_view")); 320 321 *aoout = ao; 322 PetscFunctionReturn(PETSC_SUCCESS); 323 } 324 325 /*@ 326 AOCreateMappingIS - Creates an application mapping using two index sets. 327 328 Input Parameters: 329 + isapp - index set that defines an ordering 330 - ispetsc - index set that defines another ordering, maybe NULL for identity `IS` 331 332 Output Parameter: 333 . aoout - the new application ordering 334 335 Options Database Key: 336 . -ao_view - call `AOView()` at the conclusion of `AOCreateMappingIS()` 337 338 Level: beginner 339 340 Note: 341 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. 342 Use `AOCreateBasic()` or `AOCreateBasicIS()` if they do not have holes for better performance. 343 344 .seealso: [](sec_ao), [](sec_scatter), `AOCreateBasic()`, `AOCreateMapping()`, `AODestroy()` 345 @*/ 346 PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout) 347 { 348 MPI_Comm comm; 349 const PetscInt *mypetsc, *myapp; 350 PetscInt napp, npetsc; 351 352 PetscFunctionBegin; 353 PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm)); 354 PetscCall(ISGetLocalSize(isapp, &napp)); 355 if (ispetsc) { 356 PetscCall(ISGetLocalSize(ispetsc, &npetsc)); 357 PetscCheck(napp == npetsc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local IS lengths must match"); 358 PetscCall(ISGetIndices(ispetsc, &mypetsc)); 359 } else { 360 mypetsc = NULL; 361 } 362 PetscCall(ISGetIndices(isapp, &myapp)); 363 364 PetscCall(AOCreateMapping(comm, napp, myapp, mypetsc, aoout)); 365 366 PetscCall(ISRestoreIndices(isapp, &myapp)); 367 if (ispetsc) PetscCall(ISRestoreIndices(ispetsc, &mypetsc)); 368 PetscFunctionReturn(PETSC_SUCCESS); 369 } 370