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(0); 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(0); 38 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 39 if (iascii) { 40 PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %" PetscInt_FMT "\n", aomap->N); 41 PetscViewerASCIIPrintf(viewer, " App. PETSc\n"); 42 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]]); 43 } 44 PetscFunctionReturn(0); 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(0); 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(0); 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 - Searches for the supplied application index. 126 127 Input Parameters: 128 + ao - The AOMapping 129 - index - The application index 130 131 Output Parameter: 132 . hasIndex - Flag is PETSC_TRUE if the index exists 133 134 Level: intermediate 135 136 .seealso: `AOMappingHasPetscIndex()`, `AOCreateMapping()` 137 @*/ 138 PetscErrorCode AOMappingHasApplicationIndex(AO ao, PetscInt idex, PetscBool *hasIndex) 139 { 140 AO_Mapping *aomap; 141 PetscInt *app; 142 PetscInt low, high, mid; 143 144 PetscFunctionBegin; 145 PetscValidHeaderSpecific(ao, AO_CLASSID, 1); 146 PetscValidBoolPointer(hasIndex, 3); 147 aomap = (AO_Mapping *)ao->data; 148 app = aomap->app; 149 /* Use bisection since the array is sorted */ 150 low = 0; 151 high = aomap->N - 1; 152 while (low <= high) { 153 mid = (low + high) / 2; 154 if (idex == app[mid]) break; 155 else if (idex < app[mid]) high = mid - 1; 156 else low = mid + 1; 157 } 158 if (low > high) *hasIndex = PETSC_FALSE; 159 else *hasIndex = PETSC_TRUE; 160 PetscFunctionReturn(0); 161 } 162 163 /*@ 164 AOMappingHasPetscIndex - Searches for the supplied petsc index. 165 166 Input Parameters: 167 + ao - The AOMapping 168 - index - The petsc index 169 170 Output Parameter: 171 . hasIndex - Flag is PETSC_TRUE if the index exists 172 173 Level: intermediate 174 175 .seealso: `AOMappingHasApplicationIndex()`, `AOCreateMapping()` 176 @*/ 177 PetscErrorCode AOMappingHasPetscIndex(AO ao, PetscInt idex, PetscBool *hasIndex) 178 { 179 AO_Mapping *aomap; 180 PetscInt *petsc; 181 PetscInt low, high, mid; 182 183 PetscFunctionBegin; 184 PetscValidHeaderSpecific(ao, AO_CLASSID, 1); 185 PetscValidBoolPointer(hasIndex, 3); 186 aomap = (AO_Mapping *)ao->data; 187 petsc = aomap->petsc; 188 /* Use bisection since the array is sorted */ 189 low = 0; 190 high = aomap->N - 1; 191 while (low <= high) { 192 mid = (low + high) / 2; 193 if (idex == petsc[mid]) break; 194 else if (idex < petsc[mid]) high = mid - 1; 195 else low = mid + 1; 196 } 197 if (low > high) *hasIndex = PETSC_FALSE; 198 else *hasIndex = PETSC_TRUE; 199 PetscFunctionReturn(0); 200 } 201 202 /*@C 203 AOCreateMapping - Creates a basic application mapping using two integer arrays. 204 205 Input Parameters: 206 + comm - MPI communicator that is to share AO 207 . napp - size of integer arrays 208 . myapp - integer array that defines an ordering 209 - mypetsc - integer array that defines another ordering (may be NULL to indicate the identity ordering) 210 211 Output Parameter: 212 . aoout - the new application mapping 213 214 Options Database Key: 215 . -ao_view - call AOView() at the conclusion of AOCreateMapping() 216 217 Level: beginner 218 219 Notes: 220 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. 221 Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance. 222 223 .seealso: `AOCreateBasic()`, `AOCreateBasic()`, `AOCreateMappingIS()`, `AODestroy()` 224 @*/ 225 PetscErrorCode AOCreateMapping(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout) 226 { 227 AO ao; 228 AO_Mapping *aomap; 229 PetscInt *allpetsc, *allapp; 230 PetscInt *petscPerm, *appPerm; 231 PetscInt *petsc; 232 PetscMPIInt size, rank, *lens, *disp, nnapp; 233 PetscInt N, start; 234 PetscInt i; 235 236 PetscFunctionBegin; 237 PetscValidPointer(aoout, 5); 238 *aoout = NULL; 239 PetscCall(AOInitializePackage()); 240 241 PetscCall(PetscHeaderCreate(ao, AO_CLASSID, "AO", "Application Ordering", "AO", comm, AODestroy, AOView)); 242 PetscCall(PetscNew(&aomap)); 243 PetscCall(PetscMemcpy(ao->ops, &AOps, sizeof(AOps))); 244 ao->data = (void *)aomap; 245 246 /* transmit all lengths to all processors */ 247 PetscCallMPI(MPI_Comm_size(comm, &size)); 248 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 249 PetscCall(PetscMalloc2(size, &lens, size, &disp)); 250 nnapp = napp; 251 PetscCallMPI(MPI_Allgather(&nnapp, 1, MPI_INT, lens, 1, MPI_INT, comm)); 252 N = 0; 253 for (i = 0; i < size; i++) { 254 disp[i] = N; 255 N += lens[i]; 256 } 257 aomap->N = N; 258 ao->N = N; 259 ao->n = N; 260 261 /* If mypetsc is 0 then use "natural" numbering */ 262 if (!mypetsc) { 263 start = disp[rank]; 264 PetscCall(PetscMalloc1(napp + 1, &petsc)); 265 for (i = 0; i < napp; i++) petsc[i] = start + i; 266 } else { 267 petsc = (PetscInt *)mypetsc; 268 } 269 270 /* get all indices on all processors */ 271 PetscCall(PetscMalloc4(N, &allapp, N, &appPerm, N, &allpetsc, N, &petscPerm)); 272 PetscCallMPI(MPI_Allgatherv((void *)myapp, napp, MPIU_INT, allapp, lens, disp, MPIU_INT, comm)); 273 PetscCallMPI(MPI_Allgatherv((void *)petsc, napp, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm)); 274 PetscCall(PetscFree2(lens, disp)); 275 276 /* generate a list of application and PETSc node numbers */ 277 PetscCall(PetscMalloc4(N, &aomap->app, N, &aomap->appPerm, N, &aomap->petsc, N, &aomap->petscPerm)); 278 for (i = 0; i < N; i++) { 279 appPerm[i] = i; 280 petscPerm[i] = i; 281 } 282 PetscCall(PetscSortIntWithPermutation(N, allpetsc, petscPerm)); 283 PetscCall(PetscSortIntWithPermutation(N, allapp, appPerm)); 284 /* Form sorted arrays of indices */ 285 for (i = 0; i < N; i++) { 286 aomap->app[i] = allapp[appPerm[i]]; 287 aomap->petsc[i] = allpetsc[petscPerm[i]]; 288 } 289 /* Invert petscPerm[] into aomap->petscPerm[] */ 290 for (i = 0; i < N; i++) aomap->petscPerm[petscPerm[i]] = i; 291 292 /* Form map between aomap->app[] and aomap->petsc[] */ 293 for (i = 0; i < N; i++) aomap->appPerm[i] = aomap->petscPerm[appPerm[i]]; 294 295 /* Invert appPerm[] into allapp[] */ 296 for (i = 0; i < N; i++) allapp[appPerm[i]] = i; 297 298 /* Form map between aomap->petsc[] and aomap->app[] */ 299 for (i = 0; i < N; i++) aomap->petscPerm[i] = allapp[petscPerm[i]]; 300 301 if (PetscDefined(USE_DEBUG)) { 302 /* Check that the permutations are complementary */ 303 for (i = 0; i < N; i++) PetscCheck(i == aomap->appPerm[aomap->petscPerm[i]], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid ordering"); 304 } 305 /* Cleanup */ 306 if (!mypetsc) PetscCall(PetscFree(petsc)); 307 PetscCall(PetscFree4(allapp, appPerm, allpetsc, petscPerm)); 308 309 PetscCall(AOViewFromOptions(ao, NULL, "-ao_view")); 310 311 *aoout = ao; 312 PetscFunctionReturn(0); 313 } 314 315 /*@ 316 AOCreateMappingIS - Creates a basic application ordering using two index sets. 317 318 Input Parameters: 319 + comm - MPI communicator that is to share AO 320 . isapp - index set that defines an ordering 321 - ispetsc - index set that defines another ordering, maybe NULL for identity IS 322 323 Output Parameter: 324 . aoout - the new application ordering 325 326 Options Database Key: 327 . -ao_view - call AOView() at the conclusion of AOCreateMappingIS() 328 329 Level: beginner 330 331 Notes: 332 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. 333 Use AOCreateBasic() or AOCreateBasicIS() if they do not have holes for better performance. 334 335 .seealso: `AOCreateBasic()`, `AOCreateMapping()`, `AODestroy()` 336 @*/ 337 PetscErrorCode AOCreateMappingIS(IS isapp, IS ispetsc, AO *aoout) 338 { 339 MPI_Comm comm; 340 const PetscInt *mypetsc, *myapp; 341 PetscInt napp, npetsc; 342 343 PetscFunctionBegin; 344 PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm)); 345 PetscCall(ISGetLocalSize(isapp, &napp)); 346 if (ispetsc) { 347 PetscCall(ISGetLocalSize(ispetsc, &npetsc)); 348 PetscCheck(napp == npetsc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local IS lengths must match"); 349 PetscCall(ISGetIndices(ispetsc, &mypetsc)); 350 } else { 351 mypetsc = NULL; 352 } 353 PetscCall(ISGetIndices(isapp, &myapp)); 354 355 PetscCall(AOCreateMapping(comm, napp, myapp, mypetsc, aoout)); 356 357 PetscCall(ISRestoreIndices(isapp, &myapp)); 358 if (ispetsc) PetscCall(ISRestoreIndices(ispetsc, &mypetsc)); 359 PetscFunctionReturn(0); 360 } 361