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