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

00001 
00002 #include "mpi.h"
00003 #include <cstdio>
00004 #include <iostream>
00005 #include <cassert>
00006 #include "parUtils.h"
00007 #include "handleStencils.h"
00008 
00009 int createShFnMat(double******& shFnMat) {
00010 
00011 #ifdef __USE_MG_INIT_TYPE3__
00012   createShFnMat_Type3(shFnMat);
00013 #else
00014 #ifdef __USE_MG_INIT_TYPE2__
00015   createShFnMat_Type2(shFnMat);
00016 #else
00017   createShFnMat_Type1(shFnMat);
00018 #endif
00019 #endif
00020 
00021   return 1;
00022 }
00023 
00024 int createShFnMat_Type3(double******& shFnMat) {
00025   FILE* infile;
00026   int rank;
00027   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00028 
00029   char fname[250];
00030   sprintf(fname,"ShFnStencils_%d.inp", rank);
00031   infile = fopen(fname,"r");
00032   if(!infile) {
00033     std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00034     assert(false);
00035   }
00036 
00037   typedef double* doublePtr;
00038   typedef doublePtr* double2Ptr;
00039   typedef double2Ptr* double3Ptr;
00040   typedef double3Ptr* double4Ptr;
00041   typedef double4Ptr* double5Ptr;
00042 
00043   shFnMat = new double5Ptr[8];
00044   for(unsigned int cNum = 0; cNum < 8; cNum++) {
00045     shFnMat[cNum] = new double4Ptr[18];
00046     for(unsigned int eType = 0; eType < 18; eType++) {
00047       shFnMat[cNum][eType] = new double3Ptr[8];
00048       for(unsigned int j = 0; j < 8; j++) {
00049         shFnMat[cNum][eType][j] = new double2Ptr[3];
00050         for(unsigned int m = 0; m < 3; m++) {
00051           shFnMat[cNum][eType][j][m] = new doublePtr[3];
00052           for(unsigned int n = 0; n < 3; n++) {
00053             shFnMat[cNum][eType][j][m][n] = new double[3];
00054             for(unsigned int p = 0; p < 3; p++) {
00055               fscanf(infile,"%lf",&(shFnMat[cNum][eType][j][m][n][p]));
00056             }//end p
00057           }//end n
00058         }//end m
00059       }//end j
00060     }//end etype
00061   }//end cNum
00062 
00063   fclose(infile);
00064 
00065   return 1;
00066 }
00067 
00068 int createShFnMat_Type2(double******& shFnMat) {
00069   FILE* infile;
00070   MPI_Comm comm = MPI_COMM_WORLD;
00071 
00072   int rank, npes;
00073   MPI_Comm_rank(comm, &rank);
00074   MPI_Comm_size(comm, &npes);
00075 
00076   const int THOUSAND = 1000;
00077   int numGroups = (npes/THOUSAND);
00078   if( (numGroups*THOUSAND) < npes) {
00079     numGroups++;
00080   }
00081 
00082   MPI_Comm newComm;
00083 
00084   bool* isEmptyList = new bool[npes];
00085   for(int i = 0; i < numGroups; i++) {
00086     for(int j = 0; (j < (i*THOUSAND)) && (j < npes); j++) {
00087       isEmptyList[j] = true;
00088     }
00089     for(int j = (i*THOUSAND); (j < ((i+1)*THOUSAND)) && (j < npes); j++) {
00090       isEmptyList[j] = false;
00091     }
00092     for(int j = ((i + 1)*THOUSAND); j < npes; j++) {
00093       isEmptyList[j] = true;
00094     }
00095     MPI_Comm tmpComm;
00096     par::splitComm2way(isEmptyList, &tmpComm, comm);
00097     if(!(isEmptyList[rank])) {
00098       newComm = tmpComm;
00099     }
00100   }//end for i
00101   delete [] isEmptyList;
00102 
00103 
00104   if((rank % THOUSAND) == 0) {
00105     char fname[250];
00106     sprintf(fname,"ShFnStencils_%d.inp", (rank/THOUSAND));
00107     infile = fopen(fname,"r");
00108     if(!infile) {
00109       std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00110       assert(false);
00111     }
00112   }
00113 
00114   typedef double* doublePtr;
00115   typedef doublePtr* double2Ptr;
00116   typedef double2Ptr* double3Ptr;
00117   typedef double3Ptr* double4Ptr;
00118   typedef double4Ptr* double5Ptr;
00119 
00120   shFnMat = new double5Ptr[8];
00121   for(unsigned int cNum = 0; cNum < 8; cNum++) {
00122     shFnMat[cNum] = new double4Ptr[18];
00123     for(unsigned int eType = 0; eType < 18; eType++) {
00124       shFnMat[cNum][eType] = new double3Ptr[8];
00125       for(unsigned int j = 0; j < 8; j++) {
00126         shFnMat[cNum][eType][j] = new double2Ptr[3];
00127         for(unsigned int m = 0; m < 3; m++) {
00128           shFnMat[cNum][eType][j][m] = new doublePtr[3];
00129           for(unsigned int n = 0; n < 3; n++) {
00130             shFnMat[cNum][eType][j][m][n] = new double[3];
00131             if((rank % THOUSAND) == 0) {
00132               for(unsigned int p = 0; p < 3; p++) {
00133                 fscanf(infile,"%lf",&(shFnMat[cNum][eType][j][m][n][p]));
00134               }//end p
00135             }
00136           }//end n
00137         }//end m
00138       }//end j
00139     }//end etype
00140   }//end cNum
00141 
00142   if((rank % THOUSAND) == 0) {
00143     fclose(infile);
00144   }
00145 
00146   double * tmpMat = new double[31104];
00147 
00148   if((rank % THOUSAND) == 0) {
00149     unsigned int ctr = 0;
00150     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00151       for(unsigned int eType = 0; eType < 18; eType++) {
00152         for(unsigned int j = 0; j < 8; j++) {
00153           for(unsigned int m = 0; m < 3; m++) {
00154             for(unsigned int n = 0; n < 3; n++) {
00155               for(unsigned int p = 0; p < 3; p++) {
00156                 tmpMat[ctr] = shFnMat[cNum][eType][j][m][n][p];
00157                 ctr++;
00158               }
00159             }
00160           }
00161         }
00162       }
00163     }
00164   }
00165 
00166   par::Mpi_Bcast<double>(tmpMat,31104, 0, newComm);
00167 
00168   if((rank % THOUSAND) != 0) {
00169     unsigned int ctr = 0;
00170     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00171       for(unsigned int eType = 0; eType < 18; eType++) {
00172         for(unsigned int j = 0; j < 8; j++) {
00173           for(unsigned int m = 0; m < 3; m++) {
00174             for(unsigned int n = 0; n < 3; n++) {
00175               for(unsigned int p = 0; p < 3; p++) {
00176                 shFnMat[cNum][eType][j][m][n][p] = tmpMat[ctr];
00177                 ctr++;
00178               }
00179             }
00180           }
00181         }
00182       }
00183     }
00184   }
00185 
00186   delete [] tmpMat;
00187 
00188   return 1;
00189 }//end fn.
00190 
00191 int createShFnMat_Type1(double******& shFnMat) {
00192   FILE* infile;
00193   int rank;
00194   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00195 
00196   if(!rank) {
00197     char fname[250];
00198     sprintf(fname,"ShFnStencils.inp");
00199     infile = fopen(fname,"r");
00200     if(!infile) {
00201       std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00202       assert(false);
00203     }
00204   }
00205 
00206   typedef double* doublePtr;
00207   typedef doublePtr* double2Ptr;
00208   typedef double2Ptr* double3Ptr;
00209   typedef double3Ptr* double4Ptr;
00210   typedef double4Ptr* double5Ptr;
00211 
00212   shFnMat = new double5Ptr[8];
00213   for(unsigned int cNum = 0; cNum < 8; cNum++) {
00214     shFnMat[cNum] = new double4Ptr[18];
00215     for(unsigned int eType = 0; eType < 18; eType++) {
00216       shFnMat[cNum][eType] = new double3Ptr[8];
00217       for(unsigned int j = 0; j < 8; j++) {
00218         shFnMat[cNum][eType][j] = new double2Ptr[3];
00219         for(unsigned int m = 0; m < 3; m++) {
00220           shFnMat[cNum][eType][j][m] = new doublePtr[3];
00221           for(unsigned int n = 0; n < 3; n++) {
00222             shFnMat[cNum][eType][j][m][n] = new double[3];
00223             if(!rank) {
00224               for(unsigned int p = 0; p < 3; p++) {
00225                 fscanf(infile,"%lf",&(shFnMat[cNum][eType][j][m][n][p]));
00226               }//end p
00227             }
00228           }//end n
00229         }//end m
00230       }//end j
00231     }//end etype
00232   }//end cNum
00233 
00234   if(!rank) {
00235     fclose(infile);
00236   }
00237 
00238   double * tmpMat = new double[31104];
00239 
00240   if(!rank) {
00241     unsigned int ctr = 0;
00242     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00243       for(unsigned int eType = 0; eType < 18; eType++) {
00244         for(unsigned int j = 0; j < 8; j++) {
00245           for(unsigned int m = 0; m < 3; m++) {
00246             for(unsigned int n = 0; n < 3; n++) {
00247               for(unsigned int p = 0; p < 3; p++) {
00248                 tmpMat[ctr] = shFnMat[cNum][eType][j][m][n][p];
00249                 ctr++;
00250               }
00251             }
00252           }
00253         }
00254       }
00255     }
00256   }
00257 
00258   par::Mpi_Bcast<double>(tmpMat,31104, 0, MPI_COMM_WORLD);
00259 
00260   if(rank) {
00261     unsigned int ctr = 0;
00262     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00263       for(unsigned int eType = 0; eType < 18; eType++) {
00264         for(unsigned int j = 0; j < 8; j++) {
00265           for(unsigned int m = 0; m < 3; m++) {
00266             for(unsigned int n = 0; n < 3; n++) {
00267               for(unsigned int p = 0; p < 3; p++) {
00268                 shFnMat[cNum][eType][j][m][n][p] = tmpMat[ctr];
00269                 ctr++;
00270               }
00271             }
00272           }
00273         }
00274       }
00275     }
00276   }
00277 
00278   delete [] tmpMat;
00279 
00280   return 1;
00281 }//end fn.
00282 
00283 int createGDmatType2(double ****& GDmat) {
00284 
00285 #ifdef __USE_MG_INIT_TYPE3__
00286   createGDmatType2_Type3(GDmat);
00287 #else
00288 #ifdef __USE_MG_INIT_TYPE2__
00289   createGDmatType2_Type2(GDmat);
00290 #else
00291   createGDmatType2_Type1(GDmat);
00292 #endif
00293 #endif
00294 
00295   return 1;
00296 }
00297 
00298 int createGDmatType2_Type3(double ****& GDmat) {
00299   FILE* infile;
00300   int rank;
00301   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00302 
00303   char fname[250];
00304   sprintf(fname,"GDType2Stencils_%d.inp", rank);
00305   infile = fopen(fname,"r");
00306   if(!infile) {
00307     std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00308     assert(false);
00309   }
00310 
00311   typedef double* doublePtr;
00312   typedef doublePtr* double2Ptr;
00313   typedef double2Ptr* double3Ptr;
00314 
00315   GDmat = new double3Ptr[8];
00316   for(unsigned int cNum = 0; cNum < 8; cNum++) {
00317     GDmat[cNum] = new double2Ptr[18];
00318     for(unsigned int eType = 0; eType < 18; eType++) {
00319       GDmat[cNum][eType] = new doublePtr[24];
00320       for(unsigned int i = 0; i < 24; i++) {
00321         GDmat[cNum][eType][i] = new double[24];
00322         for(unsigned int j = 0; j < 24; j++) {
00323           fscanf(infile,"%lf",&(GDmat[cNum][eType][i][j]));
00324         }
00325       }
00326     }
00327   }
00328 
00329   fclose(infile);
00330 
00331   return 1;
00332 }
00333 
00334 int createGDmatType2_Type2(double ****& GDmat) {
00335   FILE* infile;
00336   MPI_Comm comm = MPI_COMM_WORLD;
00337 
00338   int rank, npes;
00339   MPI_Comm_rank(comm, &rank);
00340   MPI_Comm_size(comm, &npes);
00341 
00342   const int THOUSAND = 1000;
00343   int numGroups = (npes/THOUSAND);
00344   if( (numGroups*THOUSAND) < npes) {
00345     numGroups++;
00346   }
00347 
00348   MPI_Comm newComm;
00349 
00350   bool* isEmptyList = new bool[npes];
00351   for(int i = 0; i < numGroups; i++) {
00352     for(int j = 0; (j < (i*THOUSAND)) && (j < npes); j++) {
00353       isEmptyList[j] = true;
00354     }
00355     for(int j = (i*THOUSAND); (j < ((i+1)*THOUSAND)) && (j < npes); j++) {
00356       isEmptyList[j] = false;
00357     }
00358     for(int j = ((i + 1)*THOUSAND); j < npes; j++) {
00359       isEmptyList[j] = true;
00360     }
00361     MPI_Comm tmpComm;
00362     par::splitComm2way(isEmptyList, &tmpComm, comm);
00363     if(!(isEmptyList[rank])) {
00364       newComm = tmpComm;
00365     }
00366   }//end for i
00367   delete [] isEmptyList;
00368 
00369 
00370   if((rank % THOUSAND) == 0) {
00371     char fname[250];
00372     sprintf(fname,"GDType2Stencils_%d.inp", (rank/THOUSAND));
00373     infile = fopen(fname,"r");
00374     if(!infile) {
00375       std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00376       assert(false);
00377     }
00378   }
00379 
00380   typedef double* doublePtr;
00381   typedef doublePtr* double2Ptr;
00382   typedef double2Ptr* double3Ptr;
00383 
00384   GDmat = new double3Ptr[8];
00385   for(unsigned int cNum = 0; cNum < 8; cNum++) {
00386     GDmat[cNum] = new double2Ptr[18];
00387     for(unsigned int eType = 0; eType < 18; eType++) {
00388       GDmat[cNum][eType] = new doublePtr[24];
00389       for(unsigned int i = 0; i < 24; i++) {
00390         GDmat[cNum][eType][i] = new double[24];
00391         if((rank % THOUSAND) == 0) {
00392           for(unsigned int j = 0; j < 24; j++) {
00393             fscanf(infile,"%lf",&(GDmat[cNum][eType][i][j]));
00394           }
00395         }
00396       }
00397     }
00398   }
00399 
00400   if((rank % THOUSAND) == 0) {
00401     fclose(infile);
00402   }
00403 
00404   double * tmpMat = new double[82944];
00405 
00406   if((rank % THOUSAND) == 0) {
00407     unsigned int ctr = 0;
00408     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00409       for(unsigned int eType = 0; eType < 18; eType++) {
00410         for(unsigned int i = 0; i < 24; i++) {
00411           for(unsigned int j = 0; j < 24; j++) {
00412             tmpMat[ctr] = GDmat[cNum][eType][i][j];
00413             ctr++;
00414           }
00415         }
00416       }
00417     }
00418   }
00419 
00420   par::Mpi_Bcast<double>(tmpMat,82944, 0, newComm);
00421 
00422   if((rank % THOUSAND) != 0) {
00423     unsigned int ctr = 0;
00424     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00425       for(unsigned int eType = 0; eType < 18; eType++) {
00426         for(unsigned int i = 0; i < 24; i++) {
00427           for(unsigned int j = 0; j < 24; j++) {
00428             GDmat[cNum][eType][i][j] = tmpMat[ctr];
00429             ctr++;
00430           }
00431         }
00432       }
00433     }
00434   }
00435 
00436   delete [] tmpMat;
00437 
00438   return 1;
00439 }//end fn.
00440 
00441 
00442 /*Type 2 Matrices: Coarse and Fine are the same.*/
00443 int createGDmatType2_Type1(double ****& GDmat) {
00444   FILE* infile;
00445   int rank;
00446   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00447 
00448   if(!rank) {
00449     char fname[250];
00450     sprintf(fname,"GDType2Stencils.inp");
00451     infile = fopen(fname,"r");
00452     if(!infile) {
00453       std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00454       assert(false);
00455     }
00456   }
00457 
00458   typedef double* doublePtr;
00459   typedef doublePtr* double2Ptr;
00460   typedef double2Ptr* double3Ptr;
00461 
00462   GDmat = new double3Ptr[8];
00463   for(unsigned int cNum = 0; cNum < 8; cNum++) {
00464     GDmat[cNum] = new double2Ptr[18];
00465     for(unsigned int eType = 0; eType < 18; eType++) {
00466       GDmat[cNum][eType] = new doublePtr[24];
00467       for(unsigned int i = 0; i < 24; i++) {
00468         GDmat[cNum][eType][i] = new double[24];
00469         if(!rank) {
00470           for(unsigned int j = 0; j < 24; j++) {
00471             fscanf(infile,"%lf",&(GDmat[cNum][eType][i][j]));
00472           }
00473         }
00474       }
00475     }
00476   }
00477 
00478   if(!rank) {
00479     fclose(infile);
00480   }
00481 
00482   double * tmpMat = new double[82944];
00483 
00484   if(!rank) {
00485     unsigned int ctr = 0;
00486     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00487       for(unsigned int eType = 0; eType < 18; eType++) {
00488         for(unsigned int i = 0; i < 24; i++) {
00489           for(unsigned int j = 0; j < 24; j++) {
00490             tmpMat[ctr] = GDmat[cNum][eType][i][j];
00491             ctr++;
00492           }
00493         }
00494       }
00495     }
00496   }
00497 
00498   par::Mpi_Bcast<double>(tmpMat,82944, 0, MPI_COMM_WORLD);
00499 
00500   if(rank) {
00501     unsigned int ctr = 0;
00502     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00503       for(unsigned int eType = 0; eType < 18; eType++) {
00504         for(unsigned int i = 0; i < 24; i++) {
00505           for(unsigned int j = 0; j < 24; j++) {
00506             GDmat[cNum][eType][i][j] = tmpMat[ctr];
00507             ctr++;
00508           }
00509         }
00510       }
00511     }
00512   }
00513 
00514   delete [] tmpMat;
00515 
00516   return 1;
00517 }//end fn.
00518 
00519 int createMmatType2(double ****& Mmat) {
00520 
00521 #ifdef __USE_MG_INIT_TYPE3__
00522   createMmatType2_Type3(Mmat);
00523 #else
00524 #ifdef __USE_MG_INIT_TYPE2__
00525   createMmatType2_Type2(Mmat);
00526 #else
00527   createMmatType2_Type1(Mmat);
00528 #endif
00529 #endif
00530 
00531   return 1;
00532 }
00533 
00534 int createMmatType2_Type3(double ****& Mmat) {
00535   FILE* infile;
00536   int rank;
00537   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00538 
00539   char fname[250];
00540   sprintf(fname,"MassType2Stencils_%d.inp", rank);
00541   infile = fopen(fname,"r");
00542   if(!infile) {
00543     std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00544     assert(false);
00545   }
00546 
00547   typedef double* doublePtr;
00548   typedef doublePtr* double2Ptr;
00549   typedef double2Ptr* double3Ptr;
00550 
00551   Mmat = new double3Ptr[8];
00552   for(unsigned int cNum = 0; cNum < 8; cNum++) {
00553     Mmat[cNum] = new double2Ptr[18];
00554     for(unsigned int eType = 0; eType < 18; eType++) {
00555       Mmat[cNum][eType] = new doublePtr[8];
00556       for(unsigned int i = 0; i < 8; i++) {
00557         Mmat[cNum][eType][i] = new double[8];
00558         for(unsigned int j = 0; j < 8; j++) {
00559           fscanf(infile,"%lf",&(Mmat[cNum][eType][i][j]));
00560         }
00561       }
00562     }
00563   }
00564 
00565   fclose(infile);
00566 
00567   return 1;
00568 }
00569 
00570 int createMmatType2_Type2(double ****& Mmat) {
00571   FILE* infile;
00572   MPI_Comm comm = MPI_COMM_WORLD;
00573 
00574   int rank, npes;
00575   MPI_Comm_rank(comm, &rank);
00576   MPI_Comm_size(comm, &npes);
00577 
00578   const int THOUSAND = 1000;
00579   int numGroups = (npes/THOUSAND);
00580   if( (numGroups*THOUSAND) < npes) {
00581     numGroups++;
00582   }
00583 
00584   MPI_Comm newComm;
00585 
00586   bool* isEmptyList = new bool[npes];
00587   for(int i = 0; i < numGroups; i++) {
00588     for(int j = 0; (j < (i*THOUSAND)) && (j < npes); j++) {
00589       isEmptyList[j] = true;
00590     }
00591     for(int j = (i*THOUSAND); (j < ((i+1)*THOUSAND)) && (j < npes); j++) {
00592       isEmptyList[j] = false;
00593     }
00594     for(int j = ((i + 1)*THOUSAND); j < npes; j++) {
00595       isEmptyList[j] = true;
00596     }
00597     MPI_Comm tmpComm;
00598     par::splitComm2way(isEmptyList, &tmpComm, comm);
00599     if(!(isEmptyList[rank])) {
00600       newComm = tmpComm;
00601     }
00602   }//end for i
00603   delete [] isEmptyList;
00604 
00605 
00606   if((rank % THOUSAND) == 0) {
00607     char fname[250];
00608     sprintf(fname,"MassType2Stencils_%d.inp", (rank/THOUSAND));
00609     infile = fopen(fname,"r");
00610     if(!infile) {
00611       std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00612       assert(false);
00613     }
00614   }
00615 
00616   typedef double* doublePtr;
00617   typedef doublePtr* double2Ptr;
00618   typedef double2Ptr* double3Ptr;
00619 
00620   Mmat = new double3Ptr[8];
00621   for(unsigned int cNum = 0; cNum < 8; cNum++) {
00622     Mmat[cNum] = new double2Ptr[18];
00623     for(unsigned int eType = 0; eType < 18; eType++) {
00624       Mmat[cNum][eType] = new doublePtr[8];
00625       for(unsigned int i = 0; i < 8; i++) {
00626         Mmat[cNum][eType][i] = new double[8];
00627         if((rank % THOUSAND) == 0) {
00628           for(unsigned int j = 0; j < 8; j++) {
00629             fscanf(infile,"%lf",&(Mmat[cNum][eType][i][j]));
00630           }
00631         }
00632       }
00633     }
00634   }
00635 
00636   if((rank % THOUSAND) == 0) {
00637     fclose(infile);
00638   }
00639 
00640   double * tmpMat = new double[9216];
00641 
00642   if((rank % THOUSAND) == 0) {
00643     unsigned int ctr = 0;
00644     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00645       for(unsigned int eType = 0; eType < 18; eType++) {
00646         for(unsigned int i = 0; i < 8; i++) {
00647           for(unsigned int j = 0; j < 8; j++) {
00648             tmpMat[ctr] = Mmat[cNum][eType][i][j];
00649             ctr++;
00650           }
00651         }
00652       }
00653     }
00654   }
00655 
00656   par::Mpi_Bcast<double>(tmpMat,9216, 0, newComm);
00657 
00658   if((rank % THOUSAND) != 0) {
00659     unsigned int ctr = 0;
00660     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00661       for(unsigned int eType = 0; eType < 18; eType++) {
00662         for(unsigned int i = 0; i < 8; i++) {
00663           for(unsigned int j = 0; j < 8; j++) {
00664             Mmat[cNum][eType][i][j] = tmpMat[ctr];
00665             ctr++;
00666           }
00667         }
00668       }
00669     }
00670   }
00671 
00672   delete [] tmpMat;
00673 
00674   return 1;
00675 }//end of function
00676 
00677 
00678 int createMmatType2_Type1(double ****& Mmat) {
00679   FILE* infile;
00680   int rank;
00681   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00682 
00683   if(!rank) {
00684     char fname[250];
00685     sprintf(fname,"MassType2Stencils.inp");
00686     infile = fopen(fname,"r");
00687     if(!infile) {
00688       std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00689       assert(false);
00690     }
00691   }
00692 
00693   typedef double* doublePtr;
00694   typedef doublePtr* double2Ptr;
00695   typedef double2Ptr* double3Ptr;
00696 
00697   Mmat = new double3Ptr[8];
00698   for(unsigned int cNum = 0; cNum < 8; cNum++) {
00699     Mmat[cNum] = new double2Ptr[18];
00700     for(unsigned int eType = 0; eType < 18; eType++) {
00701       Mmat[cNum][eType] = new doublePtr[8];
00702       for(unsigned int i = 0; i < 8; i++) {
00703         Mmat[cNum][eType][i] = new double[8];
00704         if(!rank) {
00705           for(unsigned int j = 0; j < 8; j++) {
00706             fscanf(infile,"%lf",&(Mmat[cNum][eType][i][j]));
00707           }
00708         }
00709       }
00710     }
00711   }
00712 
00713   if(!rank) {
00714     fclose(infile);
00715   }
00716 
00717   double * tmpMat = new double[9216];
00718 
00719   if(!rank) {
00720     unsigned int ctr = 0;
00721     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00722       for(unsigned int eType = 0; eType < 18; eType++) {
00723         for(unsigned int i = 0; i < 8; i++) {
00724           for(unsigned int j = 0; j < 8; j++) {
00725             tmpMat[ctr] = Mmat[cNum][eType][i][j];
00726             ctr++;
00727           }
00728         }
00729       }
00730     }
00731   }
00732 
00733   par::Mpi_Bcast<double>(tmpMat,9216, 0, MPI_COMM_WORLD);
00734 
00735   if(rank) {
00736     unsigned int ctr = 0;
00737     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00738       for(unsigned int eType = 0; eType < 18; eType++) {
00739         for(unsigned int i = 0; i < 8; i++) {
00740           for(unsigned int j = 0; j < 8; j++) {
00741             Mmat[cNum][eType][i][j] = tmpMat[ctr];
00742             ctr++;
00743           }
00744         }
00745       }
00746     }
00747   }
00748 
00749   delete [] tmpMat;
00750 
00751   return 1;
00752 }//end of function
00753 
00754 int createLmatType2(double ****& Lmat) {
00755 
00756 #ifdef __USE_MG_INIT_TYPE3__
00757   createLmatType2_Type3(Lmat);
00758 #else
00759 #ifdef __USE_MG_INIT_TYPE2__
00760   createLmatType2_Type2(Lmat);
00761 #else
00762   createLmatType2_Type1(Lmat);
00763 #endif
00764 #endif
00765 
00766   return 1;
00767 }
00768 
00769 int createLmatType2_Type3(double ****& Lmat) {
00770   FILE* infile;
00771   int rank;
00772   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00773 
00774   char fname[250];
00775   sprintf(fname,"LapType2Stencils_%d.inp", rank);
00776   infile = fopen(fname,"r");
00777   if(!infile) {
00778     std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00779     assert(false);
00780   }
00781 
00782   typedef double* doublePtr;
00783   typedef doublePtr* double2Ptr;
00784   typedef double2Ptr* double3Ptr;
00785 
00786   Lmat = new double3Ptr[8];
00787   for(unsigned int cNum = 0; cNum < 8; cNum++) {
00788     Lmat[cNum] = new double2Ptr[18];
00789     for(unsigned int eType = 0; eType < 18; eType++) {
00790       Lmat[cNum][eType] = new doublePtr[8];
00791       for(unsigned int i = 0; i < 8; i++) {
00792         Lmat[cNum][eType][i] = new double[8];
00793         for(unsigned int j = 0; j < 8; j++) {
00794           fscanf(infile,"%lf",&(Lmat[cNum][eType][i][j]));
00795         }
00796       }
00797     }
00798   }
00799 
00800   fclose(infile);
00801 
00802   return 1;
00803 }
00804 
00805 int createLmatType2_Type2(double ****& Lmat) {
00806   FILE* infile;
00807   MPI_Comm comm = MPI_COMM_WORLD;
00808 
00809   int rank, npes;
00810   MPI_Comm_rank(comm, &rank);
00811   MPI_Comm_size(comm, &npes);
00812 
00813   const int THOUSAND = 1000;
00814   int numGroups = (npes/THOUSAND);
00815   if( (numGroups*THOUSAND) < npes) {
00816     numGroups++;
00817   }
00818 
00819   MPI_Comm newComm;
00820 
00821   bool* isEmptyList = new bool[npes];
00822   for(int i = 0; i < numGroups; i++) {
00823     for(int j = 0; (j < (i*THOUSAND)) && (j < npes); j++) {
00824       isEmptyList[j] = true;
00825     }
00826     for(int j = (i*THOUSAND); (j < ((i+1)*THOUSAND)) && (j < npes); j++) {
00827       isEmptyList[j] = false;
00828     }
00829     for(int j = ((i + 1)*THOUSAND); j < npes; j++) {
00830       isEmptyList[j] = true;
00831     }
00832     MPI_Comm tmpComm;
00833     par::splitComm2way(isEmptyList, &tmpComm, comm);
00834     if(!(isEmptyList[rank])) {
00835       newComm = tmpComm;
00836     }
00837   }//end for i
00838   delete [] isEmptyList;
00839 
00840 
00841   if((rank % THOUSAND) == 0) {
00842     char fname[250];
00843     sprintf(fname,"LapType2Stencils_%d.inp", (rank/THOUSAND));
00844     infile = fopen(fname,"r");
00845     if(!infile) {
00846       std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00847       assert(false);
00848     }
00849   }
00850 
00851   typedef double* doublePtr;
00852   typedef doublePtr* double2Ptr;
00853   typedef double2Ptr* double3Ptr;
00854 
00855   Lmat = new double3Ptr[8];
00856   for(unsigned int cNum = 0; cNum < 8; cNum++) {
00857     Lmat[cNum] = new double2Ptr[18];
00858     for(unsigned int eType = 0; eType < 18; eType++) {
00859       Lmat[cNum][eType] = new doublePtr[8];
00860       for(unsigned int i = 0; i < 8; i++) {
00861         Lmat[cNum][eType][i] = new double[8];
00862         if((rank % THOUSAND) == 0) {
00863           for(unsigned int j = 0; j < 8; j++) {
00864             fscanf(infile,"%lf",&(Lmat[cNum][eType][i][j]));
00865           }
00866         }
00867       }
00868     }
00869   }
00870 
00871   if((rank % THOUSAND) == 0) {
00872     fclose(infile);
00873   }
00874 
00875   double * tmpMat = new double[9216];
00876 
00877   if((rank % THOUSAND) == 0) {
00878     unsigned int ctr = 0;
00879     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00880       for(unsigned int eType = 0; eType < 18; eType++) {
00881         for(unsigned int i = 0; i < 8; i++) {
00882           for(unsigned int j = 0; j < 8; j++) {
00883             tmpMat[ctr] = Lmat[cNum][eType][i][j];
00884             ctr++;
00885           }
00886         }
00887       }
00888     }
00889   }
00890 
00891   par::Mpi_Bcast<double>(tmpMat,9216, 0, newComm);
00892 
00893   if((rank % THOUSAND) != 0) {
00894     unsigned int ctr = 0;
00895     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00896       for(unsigned int eType = 0; eType < 18; eType++) {
00897         for(unsigned int i = 0; i < 8; i++) {
00898           for(unsigned int j = 0; j < 8; j++) {
00899             Lmat[cNum][eType][i][j] = tmpMat[ctr];
00900             ctr++;
00901           }
00902         }
00903       }
00904     }
00905   }
00906 
00907   delete [] tmpMat;
00908 
00909   return 1;
00910 }//end of function
00911 
00912 int createLmatType2_Type1(double ****& Lmat) {
00913   FILE* infile;
00914   int rank;
00915   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00916 
00917   if(!rank) {
00918     char fname[250];
00919     sprintf(fname,"LapType2Stencils.inp");
00920     infile = fopen(fname,"r");
00921     if(!infile) {
00922       std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
00923       assert(false);
00924     }
00925   }
00926 
00927   typedef double* doublePtr;
00928   typedef doublePtr* double2Ptr;
00929   typedef double2Ptr* double3Ptr;
00930 
00931   Lmat = new double3Ptr[8];
00932   for(unsigned int cNum = 0; cNum < 8; cNum++) {
00933     Lmat[cNum] = new double2Ptr[18];
00934     for(unsigned int eType = 0; eType < 18; eType++) {
00935       Lmat[cNum][eType] = new doublePtr[8];
00936       for(unsigned int i = 0; i < 8; i++) {
00937         Lmat[cNum][eType][i] = new double[8];
00938         if(!rank) {
00939           for(unsigned int j = 0; j < 8; j++) {
00940             fscanf(infile,"%lf",&(Lmat[cNum][eType][i][j]));
00941           }
00942         }
00943       }
00944     }
00945   }
00946 
00947   if(!rank) {
00948     fclose(infile);
00949   }
00950 
00951   double * tmpMat = new double[9216];
00952 
00953   if(!rank) {
00954     unsigned int ctr = 0;
00955     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00956       for(unsigned int eType = 0; eType < 18; eType++) {
00957         for(unsigned int i = 0; i < 8; i++) {
00958           for(unsigned int j = 0; j < 8; j++) {
00959             tmpMat[ctr] = Lmat[cNum][eType][i][j];
00960             ctr++;
00961           }
00962         }
00963       }
00964     }
00965   }
00966 
00967   par::Mpi_Bcast<double>(tmpMat,9216, 0, MPI_COMM_WORLD);
00968 
00969   if(rank) {
00970     unsigned int ctr = 0;
00971     for(unsigned int cNum = 0; cNum < 8; cNum++) {
00972       for(unsigned int eType = 0; eType < 18; eType++) {
00973         for(unsigned int i = 0; i < 8; i++) {
00974           for(unsigned int j = 0; j < 8; j++) {
00975             Lmat[cNum][eType][i][j] = tmpMat[ctr];
00976             ctr++;
00977           }
00978         }
00979       }
00980     }
00981   }
00982 
00983   delete [] tmpMat;
00984 
00985   return 1;
00986 }//end of function
00987 
00988 int createRHSType2(double ***& Lmat ) {
00989 
00990 #ifdef __USE_MG_INIT_TYPE3__
00991   createRHSType2_Type3(Lmat);
00992 #else
00993 #ifdef __USE_MG_INIT_TYPE2__
00994   createRHSType2_Type2(Lmat);
00995 #else
00996   createRHSType2_Type1(Lmat);
00997 #endif
00998 #endif
00999 
01000   return 1;
01001 }
01002 
01003 int createRHSType2_Type3 (double ***& Lmat ) {
01004   FILE* infile;
01005   FILE* debugfile;
01006 
01007   int rank;
01008   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
01009 
01010   char fname[250];
01011   sprintf(fname,"RHSType2Stencils_%d.inp", rank);
01012   infile = fopen(fname,"r");
01013   if(!infile) {
01014     std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
01015     assert(false);
01016   }
01017 
01018   typedef double* doublePtr;
01019   typedef doublePtr* double2Ptr;
01020 
01021   Lmat = new double2Ptr[8];
01022   for(unsigned int cNum = 0; cNum < 8; cNum++) {
01023     Lmat[cNum] = new doublePtr[18];
01024     for(unsigned int eType = 0; eType < 18; eType++) {
01025       Lmat[cNum][eType] = new double[8];
01026       for(unsigned int i = 0; i < 8; i++) {
01027         fscanf(infile,"%lf",&(Lmat[cNum][eType][i]));
01028       }
01029     }
01030   }
01031 
01032   fclose(infile);
01033 
01034   return 1;
01035 }
01036 
01037 int createRHSType2_Type2(double ***& Lmat ) {
01038   FILE* infile;
01039   FILE* debugfile;
01040   MPI_Comm comm = MPI_COMM_WORLD;
01041 
01042   int rank, npes;
01043   MPI_Comm_rank(comm, &rank);
01044   MPI_Comm_size(comm, &npes);
01045 
01046   const int THOUSAND = 1000;
01047   int numGroups = (npes/THOUSAND);
01048   if( (numGroups*THOUSAND) < npes) {
01049     numGroups++;
01050   }
01051 
01052   MPI_Comm newComm;
01053 
01054   bool* isEmptyList = new bool[npes];
01055   for(int i = 0; i < numGroups; i++) {
01056     for(int j = 0; (j < (i*THOUSAND)) && (j < npes); j++) {
01057       isEmptyList[j] = true;
01058     }
01059     for(int j = (i*THOUSAND); (j < ((i+1)*THOUSAND)) && (j < npes); j++) {
01060       isEmptyList[j] = false;
01061     }
01062     for(int j = ((i + 1)*THOUSAND); j < npes; j++) {
01063       isEmptyList[j] = true;
01064     }
01065     MPI_Comm tmpComm;
01066     par::splitComm2way(isEmptyList, &tmpComm, comm);
01067     if(!(isEmptyList[rank])) {
01068       newComm = tmpComm;
01069     }
01070   }//end for i
01071   delete [] isEmptyList;
01072 
01073 
01074   if((rank % THOUSAND) == 0) {
01075     char fname[250];
01076     sprintf(fname,"RHSType2Stencils_%d.inp", (rank/THOUSAND));
01077     infile = fopen(fname,"r");
01078     if(!infile) {
01079       std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
01080       assert(false);
01081     }
01082   }
01083 
01084   typedef double* doublePtr;
01085   typedef doublePtr* double2Ptr;
01086 
01087   Lmat = new double2Ptr[8];
01088   for(unsigned int cNum = 0; cNum < 8; cNum++) {
01089     Lmat[cNum] = new doublePtr[18];
01090     for(unsigned int eType = 0; eType < 18; eType++) {
01091       Lmat[cNum][eType] = new double[8];
01092       if((rank % THOUSAND) == 0) {
01093         for(unsigned int i = 0; i < 8; i++) {
01094           fscanf(infile,"%lf",&(Lmat[cNum][eType][i]));
01095         }
01096       }
01097     }
01098   }
01099 
01100   if((rank % THOUSAND) == 0) {
01101     fclose(infile);
01102   }
01103 
01104   const int data_size = 8*18*8;
01105   double * tmpMat = new double[data_size];
01106 
01107   if((rank % THOUSAND) == 0) {
01108     unsigned int ctr = 0;
01109     for(unsigned int cNum = 0; cNum < 8; cNum++) {
01110       for(unsigned int eType = 0; eType < 18; eType++) {
01111         for(unsigned int i = 0; i < 8; i++) {
01112           tmpMat[ctr] = Lmat[cNum][eType][i];
01113           ctr++;
01114         }
01115       }
01116     }
01117   }
01118 
01119   par::Mpi_Bcast<double>(tmpMat,data_size, 0, newComm);
01120 
01121   if((rank % THOUSAND) != 0) {
01122     unsigned int ctr = 0;
01123     for(unsigned int cNum = 0; cNum < 8; cNum++) {
01124       for(unsigned int eType = 0; eType < 18; eType++) {
01125         for(unsigned int i = 0; i < 8; i++) {
01126           Lmat[cNum][eType][i] = tmpMat[ctr];
01127           ctr++;
01128         }
01129       }
01130     }
01131   }
01132 
01133   delete [] tmpMat;
01134 
01135   return 1;
01136 }
01137 
01138 int createRHSType2_Type1 (double ***& Lmat ) {
01139   FILE* infile;
01140   FILE* debugfile;
01141 
01142   int rank;
01143   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
01144 
01145   if(!rank) {
01146     char fname[250];
01147     sprintf(fname,"RHSType2Stencils.inp");
01148     infile = fopen(fname,"r");
01149     if(!infile) {
01150       std::cout<<"The file "<<fname<<" is not good for reading."<<std::endl;
01151       assert(false);
01152     }
01153   }
01154 
01155   typedef double* doublePtr;
01156   typedef doublePtr* double2Ptr;
01157 
01158   Lmat = new double2Ptr[8];
01159   for(unsigned int cNum = 0; cNum < 8; cNum++) {
01160     Lmat[cNum] = new doublePtr[18];
01161     for(unsigned int eType = 0; eType < 18; eType++) {
01162       Lmat[cNum][eType] = new double[8];
01163       if(!rank) {
01164         for(unsigned int i = 0; i < 8; i++) {
01165           fscanf(infile,"%lf",&(Lmat[cNum][eType][i]));
01166         }
01167       }
01168     }
01169   }
01170 
01171   if(!rank) {
01172     fclose(infile);
01173   }
01174 
01175   const int data_size = 8*18*8;
01176   double * tmpMat = new double[data_size];
01177 
01178   if(!rank) {
01179     unsigned int ctr = 0;
01180     for(unsigned int cNum = 0; cNum < 8; cNum++) {
01181       for(unsigned int eType = 0; eType < 18; eType++) {
01182         for(unsigned int i = 0; i < 8; i++) {
01183           tmpMat[ctr] = Lmat[cNum][eType][i];
01184           ctr++;
01185         }
01186       }
01187     }
01188   }
01189 
01190   par::Mpi_Bcast<double>(tmpMat,data_size, 0, MPI_COMM_WORLD);
01191 
01192   if(rank) {
01193     unsigned int ctr = 0;
01194     for(unsigned int cNum = 0; cNum < 8; cNum++) {
01195       for(unsigned int eType = 0; eType < 18; eType++) {
01196         for(unsigned int i = 0; i < 8; i++) {
01197           Lmat[cNum][eType][i] = tmpMat[ctr];
01198           ctr++;
01199         }
01200       }
01201     }
01202   }
01203 
01204   delete [] tmpMat;
01205 
01206   return 1;
01207 }
01208 
01209 int destroyRHSType2 (double ***& Lmat )
01210 {  
01211   for(unsigned int cNum = 0; cNum < 8; cNum++) {
01212     for(unsigned int eType = 0; eType < 18; eType++) {
01213       delete [] (Lmat[cNum][eType]);
01214       Lmat[cNum][eType] = NULL;
01215     }
01216 
01217     delete [] (Lmat[cNum]);
01218     Lmat[cNum] = NULL;
01219   }
01220 
01221   delete [] Lmat;  
01222   Lmat = NULL;
01223 
01224   return 1;
01225 }
01226 
01227 int destroyLmatType2(double ****& Lmat ) {
01228   for(unsigned int cNum = 0; cNum < 8; cNum++) {
01229     for(unsigned int eType = 0; eType < 18; eType++) {
01230       for(unsigned int i = 0; i < 8; i++) {
01231         delete [] (Lmat[cNum][eType][i]);
01232         Lmat[cNum][eType][i] = NULL;
01233       }
01234       delete [] (Lmat[cNum][eType]);
01235       Lmat[cNum][eType] = NULL;
01236     }
01237     delete [] (Lmat[cNum]);
01238     Lmat[cNum] = NULL;
01239   }
01240 
01241   delete [] Lmat;
01242   Lmat = NULL;
01243 
01244   return 1;
01245 }//end of function
01246 
01247 int destroyMmatType2(double ****& Mmat ) {
01248   for(unsigned int cNum = 0; cNum < 8; cNum++) {
01249     for(unsigned int eType = 0; eType < 18; eType++) {
01250       for(unsigned int i = 0; i < 8; i++) {
01251         delete [] (Mmat[cNum][eType][i]);
01252         Mmat[cNum][eType][i] = NULL;
01253       }
01254       delete [] (Mmat[cNum][eType]);
01255       Mmat[cNum][eType] = NULL;
01256     }
01257     delete [] (Mmat[cNum]);
01258     Mmat[cNum] = NULL;
01259   }
01260 
01261   delete [] Mmat;
01262   Mmat = NULL;
01263 
01264   return 1;
01265 }//end of function
01266 
01267 int destroyGDmatType2(double ****& GDmat) {
01268   for(unsigned int cNum = 0; cNum < 8; cNum++) {
01269     for(unsigned int eType = 0; eType < 18; eType++) {
01270       for(unsigned int i = 0; i < 24; i++) {
01271         delete [] (GDmat[cNum][eType][i]);
01272         GDmat[cNum][eType][i] = NULL;
01273       }
01274       delete [] (GDmat[cNum][eType]);
01275       GDmat[cNum][eType] = NULL;
01276     }
01277     delete [] (GDmat[cNum]);
01278     GDmat[cNum] = NULL;
01279   }
01280 
01281   delete [] GDmat;
01282   GDmat = NULL;
01283 
01284   return 1;
01285 }//end fn.
01286 
01287 int destroyShFnMat(double******& shFnMat) {
01288   for(unsigned int cNum = 0; cNum < 8; cNum++) {
01289     for(unsigned int eType = 0; eType < 18; eType++) {
01290       for(unsigned int j = 0; j < 8; j++) {
01291         for(unsigned int m = 0; m < 3; m++) {
01292           for(unsigned int n = 0; n < 3; n++) {
01293             delete [] (shFnMat[cNum][eType][j][m][n]);
01294             shFnMat[cNum][eType][j][m][n] = NULL;
01295           }
01296           delete [] (shFnMat[cNum][eType][j][m]);
01297           shFnMat[cNum][eType][j][m] = NULL;
01298         }
01299         delete [] (shFnMat[cNum][eType][j]);
01300         shFnMat[cNum][eType][j] = NULL;
01301       }
01302       delete [] (shFnMat[cNum][eType]);
01303       shFnMat[cNum][eType] = NULL;
01304     }
01305     delete [] (shFnMat[cNum]);
01306     shFnMat[cNum] = NULL;
01307   }
01308 
01309   delete [] shFnMat;
01310   shFnMat = NULL;
01311 
01312   return 1;
01313 }//end fn.
01314 
01315 
01316