DENDRO: /net/www/csela/dendro/releases/Dendro/examples/src/backend/handleType1Stencils.C Source File

00001 
00007 #include "mpi.h"
00008 #include <cstdio>
00009 #include <cassert>
00010 #include <iostream>
00011 #include "parUtils.h"
00012 #include "handleStencils.h"
00013 
00014 /*Type 1 Matrices: Coarse is the parent of the Fine elements.*/
00015 
00016 int createMmatType1(double *****& Mmat) {
00017 
00018 #ifdef __USE_MG_INIT_TYPE3__
00019   createMmatType1_Type3(Mmat);
00020 #else
00021 #ifdef __USE_MG_INIT_TYPE2__
00022   createMmatType1_Type2(Mmat);
00023 #else
00024   createMmatType1_Type1(Mmat);
00025 #endif
00026 #endif
00027 
00028   return 1;
00029 }
00030 
00031 int createMmatType1_Type3(double *****& Mmat) {
00032   FILE* infile;
00033   int rank;
00034   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
00035 
00036   char fname[100];
00037   sprintf(fname,"MassType1Stencils_%d.inp", rank);
00038   infile = fopen(fname,"r");
00039   if(!infile) {
00040     std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00041     assert(false);
00042   }
00043 
00044   typedef double* doublePtr;
00045   typedef doublePtr* double2Ptr;
00046   typedef double2Ptr* double3Ptr;
00047   typedef double3Ptr* double4Ptr;
00048 
00049   Mmat = new double4Ptr[8];
00050   for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00051     Mmat[cNumCoarse] = new double3Ptr[18];
00052     for(unsigned int eType = 0; eType < 18; eType++) {
00053       Mmat[cNumCoarse][eType] = new double2Ptr[8];
00054       for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00055         Mmat[cNumCoarse][eType][cNumFine] = new doublePtr[8];
00056         for(unsigned int i = 0; i < 8; i++) {
00057           Mmat[cNumCoarse][eType][cNumFine][i] = new double[8];
00058           for(unsigned int j = 0; j < 8; j++) {
00059             fscanf(infile,"%lf",&(Mmat[cNumCoarse][eType][cNumFine][i][j]));
00060           }
00061         }
00062       }
00063     }
00064   }
00065 
00066   fclose(infile);
00067 
00068   return 1;
00069 }
00070 
00071 int createMmatType1_Type2(double *****& Mmat) {
00072   FILE* infile;
00073   MPI_Comm comm = MPI_COMM_WORLD;
00074 
00075   int rank, npes;
00076   MPI_Comm_rank(comm, &rank);
00077   MPI_Comm_size(comm, &npes);
00078 
00079   const int THOUSAND = 1000;
00080   int numGroups = (npes/THOUSAND);
00081   if( (numGroups*THOUSAND) < npes) {
00082     numGroups++;
00083   }
00084 
00085   MPI_Comm newComm;
00086 
00087   bool* isEmptyList = new bool[npes];
00088   for(int i = 0; i < numGroups; i++) {
00089     for(int j = 0; (j < (i*THOUSAND)) && (j < npes); j++) {
00090       isEmptyList[j] = true;
00091     }
00092     for(int j = (i*THOUSAND); (j < ((i+1)*THOUSAND)) && (j < npes); j++) {
00093       isEmptyList[j] = false;
00094     }
00095     for(int j = ((i + 1)*THOUSAND); j < npes; j++) {
00096       isEmptyList[j] = true;
00097     }
00098     MPI_Comm tmpComm;
00099     par::splitComm2way(isEmptyList, &tmpComm, comm);
00100     if(!(isEmptyList[rank])) {
00101       newComm = tmpComm;
00102     }
00103   }//end for i
00104   delete [] isEmptyList;
00105 
00106   if((rank % THOUSAND) == 0) {
00107     char fname[100];
00108     sprintf(fname,"MassType1Stencils_%d.inp",(rank/THOUSAND));
00109     infile = fopen(fname,"r");
00110     if(!infile) {
00111       std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00112       assert(false);
00113     }
00114   }
00115 
00116   typedef double* doublePtr;
00117   typedef doublePtr* double2Ptr;
00118   typedef double2Ptr* double3Ptr;
00119   typedef double3Ptr* double4Ptr;
00120 
00121   Mmat = new double4Ptr[8];
00122   for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00123     Mmat[cNumCoarse] = new double3Ptr[18];
00124     for(unsigned int eType = 0; eType < 18; eType++) {
00125       Mmat[cNumCoarse][eType] = new double2Ptr[8];
00126       for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00127         Mmat[cNumCoarse][eType][cNumFine] = new doublePtr[8];
00128         for(unsigned int i = 0; i < 8; i++) {
00129           Mmat[cNumCoarse][eType][cNumFine][i] = new double[8];
00130           if((rank % THOUSAND) == 0) {
00131             for(unsigned int j = 0; j < 8; j++) {
00132               fscanf(infile,"%lf",&(Mmat[cNumCoarse][eType][cNumFine][i][j]));
00133             }
00134           }
00135         }
00136       }
00137     }
00138   }
00139 
00140   if((rank % THOUSAND) == 0) {
00141     fclose(infile);
00142   }
00143 
00144   double * tmpMat = new double[73728];
00145 
00146   if((rank % THOUSAND) == 0) {
00147     unsigned int ctr = 0;
00148     for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00149       for(unsigned int eType = 0; eType < 18; eType++) {
00150         for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00151           for(unsigned int i = 0; i < 8; i++) {
00152             for(unsigned int j = 0; j < 8; j++) {
00153               tmpMat[ctr] = Mmat[cNumCoarse][eType][cNumFine][i][j];
00154               ctr++;
00155             }
00156           }
00157         }
00158       }
00159     }
00160   }
00161 
00162   par::Mpi_Bcast<double>(tmpMat,73728, 0, newComm);
00163 
00164   if((rank % THOUSAND) != 0) {
00165     unsigned int ctr = 0;
00166     for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00167       for(unsigned int eType = 0; eType < 18; eType++) {
00168         for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00169           for(unsigned int i = 0; i < 8; i++) {
00170             for(unsigned int j = 0; j < 8; j++) {
00171               Mmat[cNumCoarse][eType][cNumFine][i][j] = tmpMat[ctr];
00172               ctr++;
00173             }
00174           }
00175         }
00176       }
00177     }
00178   }
00179 
00180   delete [] tmpMat;
00181 
00182   return 1;
00183 }//end of function
00184 
00185 int createMmatType1_Type1(double *****& Mmat) {
00186   FILE* infile;
00187   int rank;
00188   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
00189 
00190   if(!rank) {
00191     char fname[100];
00192     sprintf(fname,"MassType1Stencils.inp");
00193     infile = fopen(fname,"r");
00194     if(!infile) {
00195       std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00196       assert(false);
00197     }
00198   }
00199 
00200   typedef double* doublePtr;
00201   typedef doublePtr* double2Ptr;
00202   typedef double2Ptr* double3Ptr;
00203   typedef double3Ptr* double4Ptr;
00204 
00205   Mmat = new double4Ptr[8];
00206   for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00207     Mmat[cNumCoarse] = new double3Ptr[18];
00208     for(unsigned int eType = 0; eType < 18; eType++) {
00209       Mmat[cNumCoarse][eType] = new double2Ptr[8];
00210       for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00211         Mmat[cNumCoarse][eType][cNumFine] = new doublePtr[8];
00212         for(unsigned int i = 0; i < 8; i++) {
00213           Mmat[cNumCoarse][eType][cNumFine][i] = new double[8];
00214           if(!rank) {
00215             for(unsigned int j = 0; j < 8; j++) {
00216               fscanf(infile,"%lf",&(Mmat[cNumCoarse][eType][cNumFine][i][j]));
00217             }
00218           }
00219         }
00220       }
00221     }
00222   }
00223 
00224   if(!rank) {
00225     fclose(infile);
00226   }
00227 
00228   double * tmpMat = new double[73728];
00229 
00230   if(!rank) {
00231     unsigned int ctr = 0;
00232     for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00233       for(unsigned int eType = 0; eType < 18; eType++) {
00234         for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00235           for(unsigned int i = 0; i < 8; i++) {
00236             for(unsigned int j = 0; j < 8; j++) {
00237               tmpMat[ctr] = Mmat[cNumCoarse][eType][cNumFine][i][j];
00238               ctr++;
00239             }
00240           }
00241         }
00242       }
00243     }
00244   }
00245 
00246   par::Mpi_Bcast<double>(tmpMat,73728, 0, MPI_COMM_WORLD);
00247 
00248   if(rank) {
00249     unsigned int ctr = 0;
00250     for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00251       for(unsigned int eType = 0; eType < 18; eType++) {
00252         for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00253           for(unsigned int i = 0; i < 8; i++) {
00254             for(unsigned int j = 0; j < 8; j++) {
00255               Mmat[cNumCoarse][eType][cNumFine][i][j] = tmpMat[ctr];
00256               ctr++;
00257             }
00258           }
00259         }
00260       }
00261     }
00262   }
00263 
00264   delete [] tmpMat;
00265 
00266   return 1;
00267 }//end of function
00268 
00269 int createLmatType1(double *****& Lmat) {
00270 
00271 #ifdef __USE_MG_INIT_TYPE3__
00272   createLmatType1_Type3(Lmat);
00273 #else
00274 #ifdef __USE_MG_INIT_TYPE2__
00275   createLmatType1_Type2(Lmat);
00276 #else
00277   createLmatType1_Type1(Lmat);
00278 #endif
00279 #endif
00280 
00281   return 1;
00282 }
00283 
00284 int createLmatType1_Type3(double *****& Lmat) {
00285   FILE* infile;
00286   int rank;
00287   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
00288 
00289   char fname[100];
00290   sprintf(fname,"LapType1Stencils_%d.inp", rank);       
00291   infile = fopen(fname,"r");
00292   if(!infile) {
00293     std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00294     assert(false);
00295   }
00296 
00297   typedef double* doublePtr;
00298   typedef doublePtr* double2Ptr;
00299   typedef double2Ptr* double3Ptr;
00300   typedef double3Ptr* double4Ptr;
00301 
00302   Lmat = new double4Ptr[8];
00303   for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00304     Lmat[cNumCoarse] = new double3Ptr[18];
00305     for(unsigned int eType = 0; eType < 18; eType++) {
00306       Lmat[cNumCoarse][eType] = new double2Ptr[8];
00307       for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00308         Lmat[cNumCoarse][eType][cNumFine] = new doublePtr[8];
00309         for(unsigned int i = 0; i < 8; i++) {
00310           Lmat[cNumCoarse][eType][cNumFine][i] = new double[8];
00311           for(unsigned int j = 0; j < 8; j++) {
00312             fscanf(infile,"%lf",&(Lmat[cNumCoarse][eType][cNumFine][i][j]));
00313           }
00314         }
00315       }
00316     }
00317   }
00318 
00319   fclose(infile);
00320 
00321   return 1;
00322 }
00323 
00324 int createLmatType1_Type2(double *****& Lmat) {
00325   FILE* infile;
00326   MPI_Comm comm = MPI_COMM_WORLD;
00327 
00328   int rank, npes;
00329   MPI_Comm_rank(comm, &rank);
00330   MPI_Comm_size(comm, &npes);
00331 
00332   const int THOUSAND = 1000;
00333   int numGroups = (npes/THOUSAND);
00334   if( (numGroups*THOUSAND) < npes) {
00335     numGroups++;
00336   }
00337 
00338   MPI_Comm newComm;
00339 
00340   bool* isEmptyList = new bool[npes];
00341   for(int i = 0; i < numGroups; i++) {
00342     for(int j = 0; (j < (i*THOUSAND)) && (j < npes); j++) {
00343       isEmptyList[j] = true;
00344     }
00345     for(int j = (i*THOUSAND); (j < ((i+1)*THOUSAND)) && (j < npes); j++) {
00346       isEmptyList[j] = false;
00347     }
00348     for(int j = ((i + 1)*THOUSAND); j < npes; j++) {
00349       isEmptyList[j] = true;
00350     }
00351     MPI_Comm tmpComm;
00352     par::splitComm2way(isEmptyList, &tmpComm, comm);
00353     if(!(isEmptyList[rank])) {
00354       newComm = tmpComm;
00355     }
00356   }//end for i
00357   delete [] isEmptyList;
00358 
00359   if((rank % THOUSAND) == 0) {
00360     char fname[100];
00361     sprintf(fname,"LapType1Stencils_%d.inp",(rank/THOUSAND));   
00362     infile = fopen(fname,"r");
00363     if(!infile) {
00364       std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00365       assert(false);
00366     }
00367   }
00368 
00369   typedef double* doublePtr;
00370   typedef doublePtr* double2Ptr;
00371   typedef double2Ptr* double3Ptr;
00372   typedef double3Ptr* double4Ptr;
00373 
00374   Lmat = new double4Ptr[8];
00375   for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00376     Lmat[cNumCoarse] = new double3Ptr[18];
00377     for(unsigned int eType = 0; eType < 18; eType++) {
00378       Lmat[cNumCoarse][eType] = new double2Ptr[8];
00379       for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00380         Lmat[cNumCoarse][eType][cNumFine] = new doublePtr[8];
00381         for(unsigned int i = 0; i < 8; i++) {
00382           Lmat[cNumCoarse][eType][cNumFine][i] = new double[8];
00383           if((rank % THOUSAND) == 0) {
00384             for(unsigned int j = 0; j < 8; j++) {
00385               fscanf(infile,"%lf",&(Lmat[cNumCoarse][eType][cNumFine][i][j]));
00386             }
00387           }
00388         }
00389       }
00390     }
00391   }
00392 
00393   if((rank % THOUSAND) == 0) {
00394     fclose(infile);
00395   }
00396 
00397   double * tmpMat = new double[73728];
00398 
00399   if((rank % THOUSAND) == 0) {
00400     unsigned int ctr = 0;
00401     for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00402       for(unsigned int eType = 0; eType < 18; eType++) {
00403         for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00404           for(unsigned int i = 0; i < 8; i++) {
00405             for(unsigned int j = 0; j < 8; j++) {
00406               tmpMat[ctr] = Lmat[cNumCoarse][eType][cNumFine][i][j];
00407               ctr++;
00408             }
00409           }
00410         }
00411       }
00412     }
00413   }
00414 
00415   par::Mpi_Bcast<double>(tmpMat,73728, 0, newComm);
00416 
00417   if((rank % THOUSAND) != 0) {
00418     unsigned int ctr = 0;
00419     for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00420       for(unsigned int eType = 0; eType < 18; eType++) {
00421         for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00422           for(unsigned int i = 0; i < 8; i++) {
00423             for(unsigned int j = 0; j < 8; j++) {
00424               Lmat[cNumCoarse][eType][cNumFine][i][j] = tmpMat[ctr];
00425               ctr++;
00426             }
00427           }
00428         }
00429       }
00430     }
00431   }
00432 
00433   delete [] tmpMat;
00434 
00435   return 1;
00436 }//end of function
00437 
00438 
00439 int createLmatType1_Type1(double *****& Lmat) {
00440   FILE* infile;
00441   int rank;
00442   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
00443 
00444   if(!rank) {
00445     char fname[100];
00446     sprintf(fname,"LapType1Stencils.inp");      
00447     infile = fopen(fname,"r");
00448     if(!infile) {
00449       std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00450       assert(false);
00451     }
00452   }
00453 
00454   typedef double* doublePtr;
00455   typedef doublePtr* double2Ptr;
00456   typedef double2Ptr* double3Ptr;
00457   typedef double3Ptr* double4Ptr;
00458 
00459   Lmat = new double4Ptr[8];
00460   for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00461     Lmat[cNumCoarse] = new double3Ptr[18];
00462     for(unsigned int eType = 0; eType < 18; eType++) {
00463       Lmat[cNumCoarse][eType] = new double2Ptr[8];
00464       for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00465         Lmat[cNumCoarse][eType][cNumFine] = new doublePtr[8];
00466         for(unsigned int i = 0; i < 8; i++) {
00467           Lmat[cNumCoarse][eType][cNumFine][i] = new double[8];
00468           if(!rank) {
00469             for(unsigned int j = 0; j < 8; j++) {
00470               fscanf(infile,"%lf",&(Lmat[cNumCoarse][eType][cNumFine][i][j]));
00471             }
00472           }
00473         }
00474       }
00475     }
00476   }
00477 
00478   if(!rank) {
00479     fclose(infile);
00480   }
00481 
00482   double * tmpMat = new double[73728];
00483 
00484   if(!rank) {
00485     unsigned int ctr = 0;
00486     for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00487       for(unsigned int eType = 0; eType < 18; eType++) {
00488         for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00489           for(unsigned int i = 0; i < 8; i++) {
00490             for(unsigned int j = 0; j < 8; j++) {
00491               tmpMat[ctr] = Lmat[cNumCoarse][eType][cNumFine][i][j];
00492               ctr++;
00493             }
00494           }
00495         }
00496       }
00497     }
00498   }
00499 
00500   par::Mpi_Bcast<double>(tmpMat,73728, 0, MPI_COMM_WORLD);
00501 
00502   if(rank) {
00503     unsigned int ctr = 0;
00504     for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00505       for(unsigned int eType = 0; eType < 18; eType++) {
00506         for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00507           for(unsigned int i = 0; i < 8; i++) {
00508             for(unsigned int j = 0; j < 8; j++) {
00509               Lmat[cNumCoarse][eType][cNumFine][i][j] = tmpMat[ctr];
00510               ctr++;
00511             }
00512           }
00513         }
00514       }
00515     }
00516   }
00517 
00518   delete [] tmpMat;
00519 
00520   return 1;
00521 }//end of function
00522 
00523 int destroyLmatType1(double *****& Lmat ) {
00524   for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00525     for(unsigned int eType = 0; eType < 18; eType++) {
00526       for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00527         for(unsigned int i = 0; i < 8; i++) {
00528           delete [] (Lmat[cNumCoarse][eType][cNumFine][i]);
00529           Lmat[cNumCoarse][eType][cNumFine][i] = NULL;
00530         }
00531         delete [] (Lmat[cNumCoarse][eType][cNumFine]);
00532         Lmat[cNumCoarse][eType][cNumFine] = NULL;
00533       }
00534       delete [] (Lmat[cNumCoarse][eType]);
00535       Lmat[cNumCoarse][eType] = NULL;
00536     }
00537     delete [] (Lmat[cNumCoarse]);
00538     Lmat[cNumCoarse] = NULL;
00539   }
00540 
00541   delete [] Lmat;
00542   Lmat = NULL;
00543 
00544   return 1;
00545 }//end of function
00546 
00547 int destroyMmatType1(double *****& Mmat ) {
00548   for(unsigned int cNumCoarse = 0; cNumCoarse < 8; cNumCoarse++) {
00549     for(unsigned int eType = 0; eType < 18; eType++) {
00550       for(unsigned int cNumFine = 0; cNumFine < 8; cNumFine++) {
00551         for(unsigned int i = 0; i < 8; i++) {
00552           delete [] (Mmat[cNumCoarse][eType][cNumFine][i]);
00553           Mmat[cNumCoarse][eType][cNumFine][i] = NULL;
00554         }
00555         delete [] (Mmat[cNumCoarse][eType][cNumFine]);
00556         Mmat[cNumCoarse][eType][cNumFine] = NULL;
00557       }
00558       delete [] (Mmat[cNumCoarse][eType]);
00559       Mmat[cNumCoarse][eType] = NULL;
00560     }
00561     delete [] (Mmat[cNumCoarse]);
00562     Mmat[cNumCoarse] = NULL;
00563   }
00564 
00565   delete [] Mmat;
00566   Mmat = NULL;
00567 
00568   return 1;
00569 }//end of function
00570 
00571 
00572