00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 static char const Id[] = "$Id: reorder.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027 #include "init.h"
00028 #include "util/util.h"
00029
00030 int CheckPartition(GridData *gdata,int parts)
00031 {
00032 MagparFunctionLogBegin;
00033
00034 int rank,size;
00035 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00036 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00037
00038 if (rank) {
00039 MagparFunctionLogReturn(0);
00040 }
00041
00042 int *nvertloc,*neleloc;
00043 ierr = PetscMalloc((parts+1)*sizeof(int),&nvertloc);CHKERRQ(ierr);
00044 ierr = PetscMalloc((parts+1)*sizeof(int),&neleloc);CHKERRQ(ierr);
00045
00046 for (int i=0;i<parts+1;i++) {
00047 nvertloc[i]=0;
00048 neleloc[i]=0;
00049 }
00050
00051
00052
00053 for (int i=0;i<gdata->n_vert;i++) {
00054 assert(gdata->vertnewproc[i]<parts);
00055 nvertloc[gdata->vertnewproc[i]+1]++;
00056 }
00057 for (int i=0;i<gdata->n_ele;i++) {
00058 assert(gdata->elenewproc[i]<parts);
00059 neleloc[gdata->elenewproc[i]+1]++;
00060 }
00061 for (int i=1;i<parts+1;i++) {
00062 nvertloc[i] += nvertloc[i-1];
00063 neleloc[i] += neleloc[i-1];
00064 }
00065
00066 PetscPrintf(PETSC_COMM_WORLD,"mesh partitioning:\n");
00067 for (int i=0;i<parts;i++) {
00068 PetscPrintf(PETSC_COMM_WORLD,
00069 "proc. %2i: %7i elements %7i vertices\n",
00070 i,
00071 neleloc[i+1]-neleloc[i],
00072 nvertloc[i+1]-nvertloc[i]
00073 );
00074 }
00075
00076 assert(nvertloc[0]==0);
00077 assert(nvertloc[parts]==gdata->n_vert);
00078 assert(neleloc[0]==0);
00079 assert(neleloc[parts]==gdata->n_ele);
00080
00081 ierr = PetscFree(nvertloc);CHKERRQ(ierr);
00082 ierr = PetscFree(neleloc);CHKERRQ(ierr);
00083
00084
00085
00086
00087
00088
00089 IS isvert;
00090 ierr = ISCreateGeneral(PETSC_COMM_SELF,gdata->n_vert,gdata->vertl2g,&isvert);CHKERRQ(ierr);
00091 ierr = ISSetPermutation(isvert);CHKERRQ(ierr);
00092
00093 Mat Atest,Btest;
00094 matcreateseqadj(gdata->n_vert,gdata->n_ele,gdata->elevert,&Atest);
00095 ierr = MatPermute(Atest,isvert,isvert,&Btest);CHKERRQ(ierr);
00096
00097 ierr = ISDestroy(isvert);CHKERRQ(ierr);
00098
00099 PetscPrintf(PETSC_COMM_WORLD,"Adjacency matrix before partitioning:\n");
00100 ierr = matviewstruct(PETSC_COMM_SELF,Atest);CHKERRQ(ierr);
00101 PetscPrintf(PETSC_COMM_WORLD,"Adjacency matrix after partitioning:\n");
00102 ierr = matviewstruct(PETSC_COMM_SELF,Btest);CHKERRQ(ierr);
00103
00104 ierr = MatDestroy(Btest);CHKERRQ(ierr);
00105 ierr = MatDestroy(Atest);CHKERRQ(ierr);
00106
00107 MagparFunctionLogReturn(0);
00108 }
00109
00110
00111
00112 int PermuteData(GridData *gdata,IS isvert,IS isele)
00113 {
00114 MagparFunctionLogBegin;
00115
00116 int rank,size;
00117 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00118 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00119
00120 if (rank) {
00121 MagparFunctionLogReturn(0);
00122 }
00123
00124 AO t_ao;
00125
00126
00127 ierr = AOCreateBasicIS(isvert,PETSC_NULL,&t_ao);CHKERRQ(ierr);
00128
00129 ierr = AOApplicationToPetscPermuteReal(t_ao,ND,gdata->vertxyz);CHKERRQ(ierr);
00130 ierr = AOApplicationToPetsc(t_ao,NV*gdata->n_ele,(PetscInt*)gdata->elevert);CHKERRQ(ierr);
00131
00132 #ifdef ADDONS
00133 if (gdata->vertcoupl!=PETSC_NULL) {
00134 ierr = AOApplicationToPetscPermuteInt(t_ao,1,gdata->vertcoupl);CHKERRQ(ierr);
00135 ierr = AOApplicationToPetsc(t_ao,gdata->n_vert,(PetscInt*)gdata->vertcoupl);CHKERRQ(ierr);
00136 }
00137 #endif
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152 PetscReal *ta_data;
00153
00154 ierr = VecGetArray(gdata->M,&ta_data);CHKERRQ(ierr);
00155 ierr = AOApplicationToPetscPermuteReal(t_ao,ND,ta_data);CHKERRQ(ierr);
00156 ierr = VecRestoreArray(gdata->M,&ta_data);CHKERRQ(ierr);
00157
00158 #ifdef ADDONS
00159 if (gdata->VH1!=PETSC_NULL) {
00160 ierr = VecGetArray(gdata->VH1,&ta_data);CHKERRQ(ierr);
00161 ierr = AOApplicationToPetscPermuteReal(t_ao,ND,ta_data);CHKERRQ(ierr);
00162 ierr = VecRestoreArray(gdata->VH1,&ta_data);CHKERRQ(ierr);
00163 }
00164 if (gdata->VH2!=PETSC_NULL) {
00165 ierr = VecGetArray(gdata->VH2,&ta_data);CHKERRQ(ierr);
00166 ierr = AOApplicationToPetscPermuteReal(t_ao,ND,ta_data);CHKERRQ(ierr);
00167 ierr = VecRestoreArray(gdata->VH2,&ta_data);CHKERRQ(ierr);
00168 }
00169
00170
00171 if (gdata->vertcoupl!=PETSC_NULL) {
00172 for (int i=0;i<gdata->n_vert;i++) {
00173 if (gdata->vertcoupl[i]<0) continue;
00174
00175 if (gdata->vertcoupl[i]>i) {
00176 int k;
00177 k=gdata->vertcoupl[i];
00178 for (int j=i+1;j<gdata->n_vert;j++) {
00179 if (gdata->vertcoupl[j]<0) continue;
00180 if (gdata->vertcoupl[j]==k) {
00181 gdata->vertcoupl[j]=i;
00182 }
00183 else if (gdata->vertcoupl[gdata->vertcoupl[j]]==i) {
00184 gdata->vertcoupl[j]=i;
00185 }
00186 else if (gdata->vertcoupl[gdata->vertcoupl[j]]==k) {
00187 gdata->vertcoupl[j]=i;
00188 }
00189 }
00190 }
00191 assert(gdata->vertcoupl[i]<i || gdata->vertcoupl[gdata->vertcoupl[i]]==i);
00192 if (gdata->vertcoupl[i]<i) {
00193 assert(gdata->vertcoupl[gdata->vertcoupl[gdata->vertcoupl[i]]]==gdata->vertcoupl[i]);
00194 }
00195 }
00196 }
00197 #endif
00198
00199 ierr = AOPetscToApplicationPermuteInt(t_ao,1,gdata->vertl2g);CHKERRQ(ierr);
00200 ascat(&gdata->vertnewproc,gdata->n_vert,1,gdata->vertnewproc,size,0);
00201 ierr = AODestroy(t_ao);CHKERRQ(ierr);
00202
00203
00204
00205 ierr = AOCreateBasicIS(isele,PETSC_NULL,&t_ao);CHKERRQ(ierr);
00206 ierr = AOApplicationToPetscPermuteInt(t_ao,NV,gdata->elevert);CHKERRQ(ierr);
00207 ierr = AOApplicationToPetscPermuteInt(t_ao,1,gdata->eleprop);CHKERRQ(ierr);
00208
00209 ierr = AOPetscToApplicationPermuteInt(t_ao,1,gdata->elel2g);CHKERRQ(ierr);
00210 ascat(&gdata->elenewproc,gdata->n_ele,1,gdata->elenewproc,size,0);
00211 ierr = AODestroy(t_ao);CHKERRQ(ierr);
00212
00213
00214
00215 for (int i=1; i<gdata->n_vert; i++)
00216 {
00217 assert(gdata->vertnewproc[i]>=gdata->vertnewproc[i-1]);
00218 assert(gdata->vertl2g[i]==i);
00219 }
00220 for (int i=1; i<gdata->n_ele; i++)
00221 {
00222 assert(gdata->elenewproc[i]>=gdata->elenewproc[i-1]);
00223 assert(gdata->elel2g[i]==i);
00224 }
00225
00226 MagparFunctionLogReturn(0);
00227 }
00228
00229
00230
00231 int OptimizeBandwidth(GridData *gdata)
00232 {
00233 PetscInt dooptimizebw;
00234 PetscTruth flg;
00235
00236 MagparFunctionLogBegin;
00237
00238 int rank,size;
00239 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00240 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00241
00242 if (rank) {
00243 MagparFunctionLogReturn(0);
00244 }
00245
00246 ierr = PetscOptionsGetInt(PETSC_NULL,"-optimizebw",&dooptimizebw,&flg);CHKERRQ(ierr);
00247 if (flg!=PETSC_TRUE) {
00248 dooptimizebw=1;
00249 PetscPrintf(PETSC_COMM_WORLD,
00250 "Option -optimizebw not found, using default value %i (enabled)!\n",
00251 dooptimizebw
00252 );
00253 }
00254 if (dooptimizebw<0) {
00255 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"-optimizebw invalid! Exiting!\n");
00256 }
00257 if (dooptimizebw==0) {
00258 PetscPrintf(PETSC_COMM_WORLD,"Matrix bandwidth optimization not performed.\n");
00259 MagparFunctionLogReturn(0);
00260 }
00261
00262 Mat Atest;
00263 IS isvert,isrow;
00264
00265 matcreateseqadj(gdata->n_vert,gdata->n_ele,gdata->elevert,&Atest);
00266
00267 ierr = MatGetOrdering(Atest,MATORDERING_RCM,&isrow,&isvert);CHKERRQ(ierr);
00268 ierr = ISDestroy(isrow);CHKERRQ(ierr);
00269 ierr = MatDestroy(Atest);CHKERRQ(ierr);
00270
00271 int iss;
00272 ierr = ISGetSize(isvert,&iss);CHKERRQ(ierr);
00273 assert(iss==gdata->n_vert);
00274
00275 #if PETSC_VERSION >= 300
00276 const PetscInt *pind;
00277 #else
00278 PetscInt *pind;
00279 #endif
00280 ierr = ISGetIndices(isvert,&pind);CHKERRQ(ierr);
00281 for (int i=0; i<gdata->n_vert; i++) {
00282 gdata->vertl2g[i]=pind[i];
00283 }
00284 ierr = ISRestoreIndices(isvert,&pind);CHKERRQ(ierr);
00285
00286 for (int i=0; i<gdata->n_ele; i++) {
00287 gdata->elel2g[i]=i;
00288 }
00289
00290
00291 for (int i=0;i<gdata->n_vert;i++) {
00292 gdata->vertnewproc[i]=distint(gdata->n_vert,size,i);
00293 }
00294
00295
00296 for (int i=0;i<gdata->n_ele;i++) {
00297
00298 gdata->elenewproc[i]=distint(gdata->n_ele,size,i);
00299
00300
00301
00302
00303
00304 }
00305
00306
00307
00308
00309 ascat(&gdata->elel2g,gdata->n_ele,1,gdata->elenewproc,size,0);
00310
00311 IS isele;
00312 ierr = ISCreateGeneral(PETSC_COMM_SELF,gdata->n_ele,gdata->elel2g,&isele);CHKERRQ(ierr);
00313 ierr = ISSetPermutation(isele);CHKERRQ(ierr);
00314
00315 ierr = CheckPartition(gdata,size);CHKERRQ(ierr);
00316 ierr = PermuteData(gdata,isvert,isele);CHKERRQ(ierr);
00317
00318 ierr = ISDestroy(isvert);CHKERRQ(ierr);
00319 ierr = ISDestroy(isele);CHKERRQ(ierr);
00320
00321
00322 MagparFunctionLogReturn(0);
00323 }
00324
00325
00326 int TrivialPartitioning(GridData *gdata)
00327 {
00328 MagparFunctionLogBegin;
00329
00330 int rank,size;
00331 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00332 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00333
00334 if (rank) {
00335 MagparFunctionLogReturn(0);
00336 }
00337
00338 #if 1
00339
00340 for (int i=0; i<gdata->n_vert; i++) {
00341 gdata->vertl2g[i]=i;
00342 }
00343 for (int i=0; i<gdata->n_ele; i++) {
00344 gdata->elel2g[i]=i;
00345 }
00346 #else
00347
00348
00349
00350 for (int i=0; i<gdata->n_vert; i++) {
00351 gdata->vertl2g[i]=gdata->n_vert-1-i;
00352 }
00353 for (int i=0; i<gdata->n_ele; i++) {
00354 gdata->elel2g[i]=gdata->n_ele-1-i;
00355 }
00356 #endif
00357
00358
00359 for (int i=0;i<gdata->n_vert;i++) {
00360 gdata->vertnewproc[i]=distint(gdata->n_vert,size,i);
00361 }
00362
00363 for (int i=0;i<gdata->n_ele;i++) {
00364
00365 gdata->elenewproc[i]=distint(gdata->n_ele,size,i);
00366
00367
00368
00369
00370
00371 }
00372
00373
00374
00375
00376 ascat(&gdata->elel2g,gdata->n_ele,1,gdata->elenewproc,size,0);
00377
00378
00379 IS isvert,isele;
00380 ierr = ISCreateGeneral(PETSC_COMM_SELF,gdata->n_vert,gdata->vertl2g,&isvert);CHKERRQ(ierr);
00381 ierr = ISSetPermutation(isvert);CHKERRQ(ierr);
00382 ierr = ISCreateGeneral(PETSC_COMM_SELF,gdata->n_ele,gdata->elel2g,&isele);CHKERRQ(ierr);
00383 ierr = ISSetPermutation(isele);CHKERRQ(ierr);
00384
00385 ierr = CheckPartition(gdata,size);CHKERRQ(ierr);
00386 #if 0
00387
00388 ierr = PermuteData(gdata,isvert,isele);CHKERRQ(ierr);
00389 #endif
00390
00391 ierr = ISDestroy(isvert);CHKERRQ(ierr);
00392 ierr = ISDestroy(isele);CHKERRQ(ierr);
00393
00394 MagparFunctionLogReturn(0);
00395 }
00396
00397
00398
00399 int Reorder(GridData *gdata)
00400 {
00401 MagparFunctionLogBegin;
00402
00403 int rank,size;
00404 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00405 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00406
00407 if (rank) {
00408 MagparFunctionLogReturn(0);
00409 }
00410
00411 ierr = PetscMalloc(gdata->n_vert*sizeof(int),&gdata->vertl2g);CHKERRQ(ierr);
00412 ierr = PetscMalloc(gdata->n_ele*sizeof(int),&gdata->elel2g);CHKERRQ(ierr);
00413 ierr = PetscMalloc(gdata->n_vert*sizeof(int),&gdata->vertnewproc);CHKERRQ(ierr);
00414 ierr = PetscMalloc(gdata->n_ele*sizeof(int),&gdata->elenewproc);CHKERRQ(ierr);
00415
00416 ierr = TrivialPartitioning(gdata);CHKERRQ(ierr);
00417
00418
00419 ierr = OptimizeBandwidth(gdata);CHKERRQ(ierr);
00420
00421
00422 ierr = DataPartitionElementsSer(gdata);CHKERRQ(ierr);
00423
00424 MagparFunctionLogReturn(0);
00425 }
00426