1 2 /* 3 The most basic AO application ordering routines. These store the 4 entire orderings on each processor. 5 */ 6 7 #include <../src/vec/is/ao/aoimpl.h> /*I "petscao.h" I*/ 8 9 typedef struct { 10 PetscInt *app; /* app[i] is the partner for the ith PETSc slot */ 11 PetscInt *petsc; /* petsc[j] is the partner for the jth app slot */ 12 } AO_Basic; 13 14 /* 15 All processors have the same data so processor 1 prints it 16 */ 17 PetscErrorCode AOView_Basic(AO ao, PetscViewer viewer) 18 { 19 PetscMPIInt rank; 20 PetscInt i; 21 AO_Basic *aobasic = (AO_Basic *)ao->data; 22 PetscBool iascii; 23 24 PetscFunctionBegin; 25 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ao), &rank)); 26 if (rank == 0) { 27 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 28 if (iascii) { 29 PetscCall(PetscViewerASCIIPrintf(viewer, "Number of elements in ordering %" PetscInt_FMT "\n", ao->N)); 30 PetscCall(PetscViewerASCIIPrintf(viewer, "PETSc->App App->PETSc\n")); 31 for (i = 0; i < ao->N; i++) PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " %3" PetscInt_FMT " %3" PetscInt_FMT " %3" PetscInt_FMT "\n", i, aobasic->app[i], i, aobasic->petsc[i])); 32 } 33 } 34 PetscCall(PetscViewerFlush(viewer)); 35 PetscFunctionReturn(0); 36 } 37 38 PetscErrorCode AODestroy_Basic(AO ao) 39 { 40 AO_Basic *aobasic = (AO_Basic *)ao->data; 41 42 PetscFunctionBegin; 43 PetscCall(PetscFree2(aobasic->app, aobasic->petsc)); 44 PetscCall(PetscFree(aobasic)); 45 PetscFunctionReturn(0); 46 } 47 48 PetscErrorCode AOBasicGetIndices_Private(AO ao, PetscInt **app, PetscInt **petsc) 49 { 50 AO_Basic *basic = (AO_Basic *)ao->data; 51 52 PetscFunctionBegin; 53 if (app) *app = basic->app; 54 if (petsc) *petsc = basic->petsc; 55 PetscFunctionReturn(0); 56 } 57 58 PetscErrorCode AOPetscToApplication_Basic(AO ao, PetscInt n, PetscInt *ia) 59 { 60 PetscInt i, N = ao->N; 61 AO_Basic *aobasic = (AO_Basic *)ao->data; 62 63 PetscFunctionBegin; 64 for (i = 0; i < n; i++) { 65 if (ia[i] >= 0 && ia[i] < N) { 66 ia[i] = aobasic->app[ia[i]]; 67 } else { 68 ia[i] = -1; 69 } 70 } 71 PetscFunctionReturn(0); 72 } 73 74 PetscErrorCode AOApplicationToPetsc_Basic(AO ao, PetscInt n, PetscInt *ia) 75 { 76 PetscInt i, N = ao->N; 77 AO_Basic *aobasic = (AO_Basic *)ao->data; 78 79 PetscFunctionBegin; 80 for (i = 0; i < n; i++) { 81 if (ia[i] >= 0 && ia[i] < N) { 82 ia[i] = aobasic->petsc[ia[i]]; 83 } else { 84 ia[i] = -1; 85 } 86 } 87 PetscFunctionReturn(0); 88 } 89 90 PetscErrorCode AOPetscToApplicationPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array) 91 { 92 AO_Basic *aobasic = (AO_Basic *)ao->data; 93 PetscInt *temp; 94 PetscInt i, j; 95 96 PetscFunctionBegin; 97 PetscCall(PetscMalloc1(ao->N * block, &temp)); 98 for (i = 0; i < ao->N; i++) { 99 for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->petsc[i] * block + j]; 100 } 101 PetscCall(PetscArraycpy(array, temp, ao->N * block)); 102 PetscCall(PetscFree(temp)); 103 PetscFunctionReturn(0); 104 } 105 106 PetscErrorCode AOApplicationToPetscPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array) 107 { 108 AO_Basic *aobasic = (AO_Basic *)ao->data; 109 PetscInt *temp; 110 PetscInt i, j; 111 112 PetscFunctionBegin; 113 PetscCall(PetscMalloc1(ao->N * block, &temp)); 114 for (i = 0; i < ao->N; i++) { 115 for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->app[i] * block + j]; 116 } 117 PetscCall(PetscArraycpy(array, temp, ao->N * block)); 118 PetscCall(PetscFree(temp)); 119 PetscFunctionReturn(0); 120 } 121 122 PetscErrorCode AOPetscToApplicationPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array) 123 { 124 AO_Basic *aobasic = (AO_Basic *)ao->data; 125 PetscReal *temp; 126 PetscInt i, j; 127 128 PetscFunctionBegin; 129 PetscCall(PetscMalloc1(ao->N * block, &temp)); 130 for (i = 0; i < ao->N; i++) { 131 for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->petsc[i] * block + j]; 132 } 133 PetscCall(PetscArraycpy(array, temp, ao->N * block)); 134 PetscCall(PetscFree(temp)); 135 PetscFunctionReturn(0); 136 } 137 138 PetscErrorCode AOApplicationToPetscPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array) 139 { 140 AO_Basic *aobasic = (AO_Basic *)ao->data; 141 PetscReal *temp; 142 PetscInt i, j; 143 144 PetscFunctionBegin; 145 PetscCall(PetscMalloc1(ao->N * block, &temp)); 146 for (i = 0; i < ao->N; i++) { 147 for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->app[i] * block + j]; 148 } 149 PetscCall(PetscArraycpy(array, temp, ao->N * block)); 150 PetscCall(PetscFree(temp)); 151 PetscFunctionReturn(0); 152 } 153 154 static struct _AOOps AOOps_Basic = { 155 PetscDesignatedInitializer(view, AOView_Basic), 156 PetscDesignatedInitializer(destroy, AODestroy_Basic), 157 PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_Basic), 158 PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_Basic), 159 PetscDesignatedInitializer(petsctoapplicationpermuteint, AOPetscToApplicationPermuteInt_Basic), 160 PetscDesignatedInitializer(applicationtopetscpermuteint, AOApplicationToPetscPermuteInt_Basic), 161 PetscDesignatedInitializer(petsctoapplicationpermutereal, AOPetscToApplicationPermuteReal_Basic), 162 PetscDesignatedInitializer(applicationtopetscpermutereal, AOApplicationToPetscPermuteReal_Basic), 163 }; 164 165 PETSC_EXTERN PetscErrorCode AOCreate_Basic(AO ao) 166 { 167 AO_Basic *aobasic; 168 PetscMPIInt size, rank, count, *lens, *disp; 169 PetscInt napp, *allpetsc, *allapp, ip, ia, N, i, *petsc = NULL, start; 170 IS isapp = ao->isapp, ispetsc = ao->ispetsc; 171 MPI_Comm comm; 172 const PetscInt *myapp, *mypetsc = NULL; 173 174 PetscFunctionBegin; 175 /* create special struct aobasic */ 176 PetscCall(PetscNew(&aobasic)); 177 ao->data = (void *)aobasic; 178 PetscCall(PetscMemcpy(ao->ops, &AOOps_Basic, sizeof(struct _AOOps))); 179 PetscCall(PetscObjectChangeTypeName((PetscObject)ao, AOBASIC)); 180 181 PetscCall(ISGetLocalSize(isapp, &napp)); 182 PetscCall(ISGetIndices(isapp, &myapp)); 183 184 PetscCall(PetscMPIIntCast(napp, &count)); 185 186 /* transmit all lengths to all processors */ 187 PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm)); 188 PetscCallMPI(MPI_Comm_size(comm, &size)); 189 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 190 PetscCall(PetscMalloc2(size, &lens, size, &disp)); 191 PetscCallMPI(MPI_Allgather(&count, 1, MPI_INT, lens, 1, MPI_INT, comm)); 192 N = 0; 193 for (i = 0; i < size; i++) { 194 PetscCall(PetscMPIIntCast(N, disp + i)); /* = sum(lens[j]), j< i */ 195 N += lens[i]; 196 } 197 ao->N = N; 198 ao->n = N; 199 200 /* If mypetsc is 0 then use "natural" numbering */ 201 if (napp) { 202 if (!ispetsc) { 203 start = disp[rank]; 204 PetscCall(PetscMalloc1(napp + 1, &petsc)); 205 for (i = 0; i < napp; i++) petsc[i] = start + i; 206 } else { 207 PetscCall(ISGetIndices(ispetsc, &mypetsc)); 208 petsc = (PetscInt *)mypetsc; 209 } 210 } 211 212 /* get all indices on all processors */ 213 PetscCall(PetscMalloc2(N, &allpetsc, N, &allapp)); 214 PetscCallMPI(MPI_Allgatherv(petsc, count, MPIU_INT, allpetsc, lens, disp, MPIU_INT, comm)); 215 PetscCallMPI(MPI_Allgatherv((void *)myapp, count, MPIU_INT, allapp, lens, disp, MPIU_INT, comm)); 216 PetscCall(PetscFree2(lens, disp)); 217 218 if (PetscDefined(USE_DEBUG)) { 219 PetscInt *sorted; 220 PetscCall(PetscMalloc1(N, &sorted)); 221 222 PetscCall(PetscArraycpy(sorted, allpetsc, N)); 223 PetscCall(PetscSortInt(N, sorted)); 224 for (i = 0; i < N; i++) PetscCheck(sorted[i] == i, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PETSc ordering requires a permutation of numbers 0 to N-1\n it is missing %" PetscInt_FMT " has %" PetscInt_FMT, i, sorted[i]); 225 226 PetscCall(PetscArraycpy(sorted, allapp, N)); 227 PetscCall(PetscSortInt(N, sorted)); 228 for (i = 0; i < N; i++) PetscCheck(sorted[i] == i, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Application ordering requires a permutation of numbers 0 to N-1\n it is missing %" PetscInt_FMT " has %" PetscInt_FMT, i, sorted[i]); 229 230 PetscCall(PetscFree(sorted)); 231 } 232 233 /* generate a list of application and PETSc node numbers */ 234 PetscCall(PetscCalloc2(N, &aobasic->app, N, &aobasic->petsc)); 235 for (i = 0; i < N; i++) { 236 ip = allpetsc[i]; 237 ia = allapp[i]; 238 /* check there are no duplicates */ 239 PetscCheck(!aobasic->app[ip], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Duplicate in PETSc ordering at position %" PetscInt_FMT ". Already mapped to %" PetscInt_FMT ", not %" PetscInt_FMT ".", i, aobasic->app[ip] - 1, ia); 240 aobasic->app[ip] = ia + 1; 241 PetscCheck(!aobasic->petsc[ia], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Duplicate in Application ordering at position %" PetscInt_FMT ". Already mapped to %" PetscInt_FMT ", not %" PetscInt_FMT ".", i, aobasic->petsc[ia] - 1, ip); 242 aobasic->petsc[ia] = ip + 1; 243 } 244 if (napp && !mypetsc) PetscCall(PetscFree(petsc)); 245 PetscCall(PetscFree2(allpetsc, allapp)); 246 /* shift indices down by one */ 247 for (i = 0; i < N; i++) { 248 aobasic->app[i]--; 249 aobasic->petsc[i]--; 250 } 251 252 PetscCall(ISRestoreIndices(isapp, &myapp)); 253 if (napp) { 254 if (ispetsc) { 255 PetscCall(ISRestoreIndices(ispetsc, &mypetsc)); 256 } else { 257 PetscCall(PetscFree(petsc)); 258 } 259 } 260 PetscFunctionReturn(0); 261 } 262 263 /*@C 264 AOCreateBasic - Creates a basic application ordering using two integer arrays. 265 266 Collective 267 268 Input Parameters: 269 + comm - MPI communicator that is to share AO 270 . napp - size of integer arrays 271 . myapp - integer array that defines an ordering 272 - mypetsc - integer array that defines another ordering (may be NULL to 273 indicate the natural ordering, that is 0,1,2,3,...) 274 275 Output Parameter: 276 . aoout - the new application ordering 277 278 Level: beginner 279 280 Notes: 281 the arrays myapp and mypetsc must contain the all the integers 0 to napp-1 with no duplicates; that is there cannot be any "holes" 282 in the indices. Use AOCreateMapping() or AOCreateMappingIS() if you wish to have "holes" in the indices. 283 284 .seealso: `AOCreateBasicIS()`, `AODestroy()`, `AOPetscToApplication()`, `AOApplicationToPetsc()` 285 @*/ 286 PetscErrorCode AOCreateBasic(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout) 287 { 288 IS isapp, ispetsc; 289 const PetscInt *app = myapp, *petsc = mypetsc; 290 291 PetscFunctionBegin; 292 PetscCall(ISCreateGeneral(comm, napp, app, PETSC_USE_POINTER, &isapp)); 293 if (mypetsc) { 294 PetscCall(ISCreateGeneral(comm, napp, petsc, PETSC_USE_POINTER, &ispetsc)); 295 } else { 296 ispetsc = NULL; 297 } 298 PetscCall(AOCreateBasicIS(isapp, ispetsc, aoout)); 299 PetscCall(ISDestroy(&isapp)); 300 if (mypetsc) PetscCall(ISDestroy(&ispetsc)); 301 PetscFunctionReturn(0); 302 } 303 304 /*@C 305 AOCreateBasicIS - Creates a basic application ordering using two index sets. 306 307 Collective on IS 308 309 Input Parameters: 310 + isapp - index set that defines an ordering 311 - ispetsc - index set that defines another ordering (may be NULL to use the 312 natural ordering) 313 314 Output Parameter: 315 . aoout - the new application ordering 316 317 Level: beginner 318 319 Notes: 320 the index sets isapp and ispetsc must contain the all the integers 0 to napp-1 (where napp is the length of the index sets) with no duplicates; 321 that is there cannot be any "holes" 322 323 .seealso: `AOCreateBasic()`, `AODestroy()` 324 @*/ 325 PetscErrorCode AOCreateBasicIS(IS isapp, IS ispetsc, AO *aoout) 326 { 327 MPI_Comm comm; 328 AO ao; 329 330 PetscFunctionBegin; 331 PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm)); 332 PetscCall(AOCreate(comm, &ao)); 333 PetscCall(AOSetIS(ao, isapp, ispetsc)); 334 PetscCall(AOSetType(ao, AOBASIC)); 335 PetscCall(AOViewFromOptions(ao, NULL, "-ao_view")); 336 *aoout = ao; 337 PetscFunctionReturn(0); 338 } 339