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