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