fmm3d_mpi: main

00023 {
00024   PetscInitialize(&argc,&argv,"options_0",NULL);  //PetscTruth flg = PETSC_FALSE;
00025 
00026   MPI_Comm comm; comm = PETSC_COMM_WORLD;
00027   int mpirank;  MPI_Comm_rank(MPI_COMM_WORLD, &mpirank);
00028   int mpisize;  MPI_Comm_size(MPI_COMM_WORLD, &mpisize);
00029   int dim = 3;
00030   srand48( (long)( time(NULL)+mpirank ) );
00031   
00032   PetscTruth flg = PETSC_FALSE;
00033   pC( PetscPrintf(MPI_COMM_WORLD, "mpisize %d\n", mpisize) );
00034   
00035   //1. random src and trg data
00036   int numsrc;  pC( PetscOptionsGetInt(PETSC_NULL, "-numsrc", &numsrc, &flg) ); pA(flg==PETSC_TRUE);
00037   int numtrg;  pC( PetscOptionsGetInt(PETSC_NULL, "-numtrg", &numtrg, &flg) ); pA(flg==PETSC_TRUE);
00038   int kt;  pC( PetscOptionsGetInt(PETSC_NULL, "-kt", &kt, &flg) ); pA(flg==PETSC_TRUE);
00039   
00040   vector<double> tmp(2);         tmp[0] = 1;     tmp[1] = 0.25; //coefs in the kernel, work for all examples
00041   Kernel3d_MPI knl(kt, tmp);
00042   
00043   int lclnumsrc = (numsrc+mpirank)/mpisize;
00044   int lclnumtrg = (numtrg+mpirank)/mpisize;
00045   
00046   //pick center
00047   double R = 0.85;
00048   double r = 1-R;
00049   double lclctr[3];
00050   for(int d=0; d<dim; d++)       lclctr[d] = R * (2.0*drand48()-1.0);
00051   
00052   //generate src and trg data
00053   Vec srcPos;  pC( VecCreateMPI(comm, lclnumsrc*dim, PETSC_DETERMINE, &srcPos) );
00054   double* srcPosarr; pC( VecGetArray(srcPos, &srcPosarr) );
00055   for(int k=0; k<lclnumsrc; k++)
00056          for(int d=0; d<dim; d++)
00057                 srcPosarr[d+dim*k] = lclctr[d] + r*(2.0*drand48()-1.0);
00058   pC( VecRestoreArray(srcPos, &srcPosarr) );
00059   
00060   Vec trgPos;  pC( VecCreateMPI(comm, lclnumtrg*dim, PETSC_DETERMINE, &trgPos) );
00061   double* trgPosarr; pC( VecGetArray(trgPos, &trgPosarr) );
00062   for(int k=0; k<lclnumtrg; k++)
00063          for(int d=0; d<dim; d++)
00064                 trgPosarr[d+dim*k] = lclctr[d] + r*(2.0*drand48()-1.0);
00065   pC( VecRestoreArray(trgPos, &trgPosarr) );
00066   
00067   int srcDOF = knl.srcDOF();
00068   int trgDOF = knl.trgDOF();
00069   Vec srcDen;  pC( VecCreateMPI(comm, lclnumsrc*srcDOF, PETSC_DETERMINE, &srcDen) );
00070   double* srcDenarr; pC( VecGetArray(srcDen, &srcDenarr) );
00071   for(int k=0; k<lclnumsrc*srcDOF; k++)
00072          srcDenarr[k] = drand48();
00073   pC( VecRestoreArray(srcDen, &srcDenarr) );
00074   
00075   Vec trgVal;  pC( VecCreateMPI(comm, lclnumtrg*trgDOF, PETSC_DETERMINE, &trgVal) );
00076   double zero = 0;  pC( VecSet(trgVal, zero) );
00077   
00078   //2. allocate fmm 
00079   clock_t ck0, ck1;
00080   
00081   FMM3d_MPI* fmm = new FMM3d_MPI("fmm3d_MPI_");
00082   fmm->srcPos()=srcPos;  fmm->srcNor()=srcPos;
00083   fmm->trgPos()=trgPos;
00084   fmm->ctr() = Point3(0,0,0); // CENTER OF THE TOPLEVEL BOX
00085   fmm->rootLevel() = 0;         // 2^(-rootlvl) is the RADIUS OF THE TOPLEVEL BOX
00086   fmm->knl() = knl;
00087   //setup
00088   ck0 = clock();
00089   pC( fmm->setup() );
00090   ck1 = clock();  
00091   pC( PetscPrintf(MPI_COMM_WORLD, "fmm setup used %e secs\n", double(ck1-ck0)/CLOCKS_PER_SEC) );
00092   
00093   //3. run fmm
00094   for(int i=0; i<3; i++) {
00095          ck0 = clock();
00096          pC( fmm->evaluate(srcDen, trgVal) );
00097          //pC( MPI_Barrier(MPI_COMM_WORLD) );  //check Vincnt, Votcnt
00098          ck1 = clock();
00099          pC( PetscPrintf(MPI_COMM_WORLD, "fmm evaluate used %e secs\n", double(ck1-ck0)/CLOCKS_PER_SEC) );
00100   }
00101   
00102   //4. check
00103   ck0 = clock();
00104   double rerr;
00105   pC( fmm->check(srcDen, trgVal, 20, rerr) );
00106   pC( PetscPrintf(MPI_COMM_WORLD, "relative %e\n", rerr) );
00107   ck1 = clock();
00108   pC( PetscPrintf(MPI_COMM_WORLD, "fmm check used %e secs\n", double(ck1-ck0)/CLOCKS_PER_SEC) );
00109   //extra
00110   //{ PetscViewer viewer;  pC( PetscViewerASCIIOpen(comm, "srcPos.out", &viewer) );  pC( VecView(srcPos, viewer) );  pC( PetscViewerDestroy(viewer) ); }
00111   //{ PetscViewer viewer;  pC( PetscViewerASCIIOpen(comm, "trgPos.out", &viewer) );  pC( VecView(trgPos, viewer) );  pC( PetscViewerDestroy(viewer) ); }
00112   
00113   delete fmm;
00114   
00115   pC( VecDestroy(srcPos) );
00116   pC( VecDestroy(trgPos) );
00117   pC( VecDestroy(srcDen) );
00118   pC( VecDestroy(trgVal) );
00119   
00120   PetscFinalize();
00121   return 0;
00122 }