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