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 static 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(PETSC_SUCCESS); 36 } 37 38 static 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(PETSC_SUCCESS); 46 } 47 48 static PetscErrorCode AOPetscToApplication_Basic(AO ao, PetscInt n, PetscInt *ia) 49 { 50 PetscInt i, N = ao->N; 51 AO_Basic *aobasic = (AO_Basic *)ao->data; 52 53 PetscFunctionBegin; 54 for (i = 0; i < n; i++) { 55 if (ia[i] >= 0 && ia[i] < N) { 56 ia[i] = aobasic->app[ia[i]]; 57 } else { 58 ia[i] = -1; 59 } 60 } 61 PetscFunctionReturn(PETSC_SUCCESS); 62 } 63 64 static PetscErrorCode AOApplicationToPetsc_Basic(AO ao, PetscInt n, PetscInt *ia) 65 { 66 PetscInt i, N = ao->N; 67 AO_Basic *aobasic = (AO_Basic *)ao->data; 68 69 PetscFunctionBegin; 70 for (i = 0; i < n; i++) { 71 if (ia[i] >= 0 && ia[i] < N) { 72 ia[i] = aobasic->petsc[ia[i]]; 73 } else { 74 ia[i] = -1; 75 } 76 } 77 PetscFunctionReturn(PETSC_SUCCESS); 78 } 79 80 static PetscErrorCode AOPetscToApplicationPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array) 81 { 82 AO_Basic *aobasic = (AO_Basic *)ao->data; 83 PetscInt *temp; 84 PetscInt i, j; 85 86 PetscFunctionBegin; 87 PetscCall(PetscMalloc1(ao->N * block, &temp)); 88 for (i = 0; i < ao->N; i++) { 89 for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->petsc[i] * block + j]; 90 } 91 PetscCall(PetscArraycpy(array, temp, ao->N * block)); 92 PetscCall(PetscFree(temp)); 93 PetscFunctionReturn(PETSC_SUCCESS); 94 } 95 96 static PetscErrorCode AOApplicationToPetscPermuteInt_Basic(AO ao, PetscInt block, PetscInt *array) 97 { 98 AO_Basic *aobasic = (AO_Basic *)ao->data; 99 PetscInt *temp; 100 PetscInt i, j; 101 102 PetscFunctionBegin; 103 PetscCall(PetscMalloc1(ao->N * block, &temp)); 104 for (i = 0; i < ao->N; i++) { 105 for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->app[i] * block + j]; 106 } 107 PetscCall(PetscArraycpy(array, temp, ao->N * block)); 108 PetscCall(PetscFree(temp)); 109 PetscFunctionReturn(PETSC_SUCCESS); 110 } 111 112 static PetscErrorCode AOPetscToApplicationPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array) 113 { 114 AO_Basic *aobasic = (AO_Basic *)ao->data; 115 PetscReal *temp; 116 PetscInt i, j; 117 118 PetscFunctionBegin; 119 PetscCall(PetscMalloc1(ao->N * block, &temp)); 120 for (i = 0; i < ao->N; i++) { 121 for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->petsc[i] * block + j]; 122 } 123 PetscCall(PetscArraycpy(array, temp, ao->N * block)); 124 PetscCall(PetscFree(temp)); 125 PetscFunctionReturn(PETSC_SUCCESS); 126 } 127 128 static PetscErrorCode AOApplicationToPetscPermuteReal_Basic(AO ao, PetscInt block, PetscReal *array) 129 { 130 AO_Basic *aobasic = (AO_Basic *)ao->data; 131 PetscReal *temp; 132 PetscInt i, j; 133 134 PetscFunctionBegin; 135 PetscCall(PetscMalloc1(ao->N * block, &temp)); 136 for (i = 0; i < ao->N; i++) { 137 for (j = 0; j < block; j++) temp[i * block + j] = array[aobasic->app[i] * block + j]; 138 } 139 PetscCall(PetscArraycpy(array, temp, ao->N * block)); 140 PetscCall(PetscFree(temp)); 141 PetscFunctionReturn(PETSC_SUCCESS); 142 } 143 144 static const struct _AOOps AOOps_Basic = { 145 PetscDesignatedInitializer(view, AOView_Basic), 146 PetscDesignatedInitializer(destroy, AODestroy_Basic), 147 PetscDesignatedInitializer(petsctoapplication, AOPetscToApplication_Basic), 148 PetscDesignatedInitializer(applicationtopetsc, AOApplicationToPetsc_Basic), 149 PetscDesignatedInitializer(petsctoapplicationpermuteint, AOPetscToApplicationPermuteInt_Basic), 150 PetscDesignatedInitializer(applicationtopetscpermuteint, AOApplicationToPetscPermuteInt_Basic), 151 PetscDesignatedInitializer(petsctoapplicationpermutereal, AOPetscToApplicationPermuteReal_Basic), 152 PetscDesignatedInitializer(applicationtopetscpermutereal, AOApplicationToPetscPermuteReal_Basic), 153 }; 154 155 PETSC_INTERN PetscErrorCode AOCreate_Basic(AO ao) 156 { 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(PetscNew(&aobasic)); 167 ao->data = (void *)aobasic; 168 ao->ops[0] = AOOps_Basic; 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 for (i = 0; i < N; i++) { 226 ip = allpetsc[i]; 227 ia = allapp[i]; 228 /* check there are no duplicates */ 229 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); 230 aobasic->app[ip] = ia + 1; 231 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); 232 aobasic->petsc[ia] = ip + 1; 233 } 234 if (napp && !mypetsc) PetscCall(PetscFree(petsc)); 235 PetscCall(PetscFree2(allpetsc, allapp)); 236 /* shift indices down by one */ 237 for (i = 0; i < N; i++) { 238 aobasic->app[i]--; 239 aobasic->petsc[i]--; 240 } 241 242 PetscCall(ISRestoreIndices(isapp, &myapp)); 243 if (napp) { 244 if (ispetsc) { 245 PetscCall(ISRestoreIndices(ispetsc, &mypetsc)); 246 } else { 247 PetscCall(PetscFree(petsc)); 248 } 249 } 250 PetscFunctionReturn(PETSC_SUCCESS); 251 } 252 253 /*@C 254 AOCreateBasic - Creates a basic application ordering using two integer arrays. 255 256 Collective 257 258 Input Parameters: 259 + comm - MPI communicator that is to share `AO` 260 . napp - size of integer arrays 261 . myapp - integer array that defines an ordering 262 - mypetsc - integer array that defines another ordering (may be NULL to 263 indicate the natural ordering, that is 0,1,2,3,...) 264 265 Output Parameter: 266 . aoout - the new application ordering 267 268 Level: beginner 269 270 Note: 271 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" 272 in the indices. Use `AOCreateMapping()` or `AOCreateMappingIS()` if you wish to have "holes" in the indices. 273 274 .seealso: [](sec_ao), [](sec_scatter), `AO`, `AOCreateBasicIS()`, `AODestroy()`, `AOPetscToApplication()`, `AOApplicationToPetsc()` 275 @*/ 276 PetscErrorCode AOCreateBasic(MPI_Comm comm, PetscInt napp, const PetscInt myapp[], const PetscInt mypetsc[], AO *aoout) 277 { 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(PETSC_SUCCESS); 292 } 293 294 /*@C 295 AOCreateBasicIS - Creates a basic application ordering using two `IS` index sets. 296 297 Collective 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 Note: 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: [](sec_ao), [](sec_scatter), `IS`, `AO`, `AOCreateBasic()`, `AODestroy()` 314 @*/ 315 PetscErrorCode AOCreateBasicIS(IS isapp, IS ispetsc, AO *aoout) 316 { 317 MPI_Comm comm; 318 AO ao; 319 320 PetscFunctionBegin; 321 PetscCall(PetscObjectGetComm((PetscObject)isapp, &comm)); 322 PetscCall(AOCreate(comm, &ao)); 323 PetscCall(AOSetIS(ao, isapp, ispetsc)); 324 PetscCall(AOSetType(ao, AOBASIC)); 325 PetscCall(AOViewFromOptions(ao, NULL, "-ao_view")); 326 *aoout = ao; 327 PetscFunctionReturn(PETSC_SUCCESS); 328 } 329