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