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