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