Main Page | Class Hierarchy | Class List | File List | Class Members | File Members

mdp_lattice.h

Go to the documentation of this file.
00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 #define MDP_LATTICE
00014 
00015 const long NOWHERE=INT_MAX;
00016                    
00017 class mdp_site;
00018 
00031 class mdp_lattice {
00032  public:
00033   int ndim;        /* number of dimensions       */
00034   int ndir;        /* number of directions       */
00035   int next_next;   /* 1 or 2 is the thikness of the boudary */
00036   int *nx;         /* box containing the lattice */
00037   long nvol;       /* local volume               */
00038   long nvol_gl;    /* global volume              */
00039   long nvol_in;    /* internal volume            */
00040   long *gl;        /* map local to global        */
00041   long *lg;        /* map global to local        */
00042   FILE *lg_file;   /* temporary file to store lg if not enough memory */
00043   long **up;       /* move up in local index     */
00044   long **dw;       /* move dw in local index     */
00045   int  **co;       /* coordinate x in local idx  */
00046   int  *wh;        /* in which process? is idx   */
00047   int  *parity;    /* parity of local mdp_site idx   */
00048   long start[_NprocMax_][2]; 
00049   long stop[_NprocMax_][2];
00050   long len_to_send[_NprocMax_][2];
00051   long *to_send[_NprocMax_];
00052   bool local_random_generator;
00053  private:
00054   mdp_prng *random_obj;
00055   long random_seed;
00056   // ////////////////////////////////////////////////////
00057   // share information with other processors
00058   // ////////////////////////////////////////////////////
00059   void communicate_results_to_all_processes() {
00060     long buffer[2];
00061     long *dynamic_buffer;
00062     long length;
00063     int dp, process, np, idx;
00064     mdp_request request;
00065 
00066     // sending length ///////////////////////////
00067     for(dp=1; dp<Nproc; dp++) {
00068       process=(ME+dp) % Nproc;
00069       for(np=0; np<2; np++)
00070         buffer[np]=stop[process][np]-start[process][np];
00071       mpi.put(buffer, 2, process, request);
00072       process=(ME-dp+Nproc) % Nproc;
00073       mpi.get(len_to_send[process], 2, process);
00074       mpi.wait(request);
00075       process=(ME+dp) % Nproc;
00076       length=stop[process][1]-start[process][0];
00077       dynamic_buffer= new long[length];
00078       for(idx=0; idx<length; idx++)
00079         dynamic_buffer[idx]=gl[start[process][0]+idx];
00080       mpi.put(dynamic_buffer, length, process, request);
00081       process=(ME-dp+Nproc) % Nproc;
00082       length=len_to_send[process][0]+len_to_send[process][1];
00083       to_send[process]=new long[length];
00084       mpi.get(to_send[process], length, process);
00085       for(idx=0; idx<length; idx++) 
00086         to_send[process][idx]=local(to_send[process][idx]);
00087       mpi.wait(request);
00088       delete[] dynamic_buffer;
00089     }
00090   }
00091  public:
00092   int (*where)(int*, int, int*);
00093   void (*neighbour)(int, int*, int*, int*, int, int*);
00094   inline long global_coordinate(int *x) {
00095     long  global_idx=0;
00096     int   mu;
00097     for(mu=0; mu<ndim-1; mu++) global_idx=(global_idx+x[mu])*nx[mu+1];
00098     return global_idx+x[ndim-1];
00099   }
00100   inline void global_coordinate(long global_idx, int *x) {
00101     int mu;
00102     for(mu=ndim-1; mu>0; mu--) {
00103       x[mu]=global_idx % nx[mu];
00104       global_idx=(global_idx-x[mu])/nx[mu];
00105     }
00106     x[0]=global_idx;
00107   }
00108   inline int compute_parity(int *x) {
00109     int mu=0;
00110     int p=0;
00111     for(mu=0; mu<ndim; mu++) p=p+x[mu];
00112     return (p % 2);
00113   }
00114   mdp_lattice() {
00115     nvol=0;
00116     random_obj=0;
00117   }
00126   mdp_lattice(int ndim_, 
00127               int nx_[],
00128               int (*where_)(int*, int, int*)=default_partitioning0, 
00129               void (*neighbour_)(int,int*,int*,int*,int,int*) 
00130               =torus_topology,
00131               long random_seed_=0,
00132               int next_next_=1,
00133               bool local_random_=true) {
00134     nvol=0;
00135     random_obj=0;
00136     allocate_lattice(ndim_,ndim_,nx_,where_,neighbour_,
00137                      random_seed_, next_next_, local_random_);
00138   }    
00140   mdp_lattice(int ndim_, 
00141               int ndir_,
00142               int nx_[],
00143               int (*where_)(int*, int, int*)=default_partitioning0,
00144               void (*neighbour_)(int,int*,int*,int*,int,int*) 
00145               =torus_topology,
00146               long random_seed_=0,
00147               int next_next_=1,
00148               bool local_random_=true) {
00149     nvol=0;
00150     random_obj=0;
00151     allocate_lattice(ndim_,ndir_,nx_,where_,neighbour_,
00152                      random_seed_, next_next_, local_random_);
00153   }   
00162   void allocate_lattice(int ndim_, 
00163                         int nx_[],
00164                         int (*where_)(int*, int, int*)=default_partitioning0,
00165                         void (*neighbour_)(int,int*,int*,int*,int,int*) 
00166                         =torus_topology,
00167                         long random_seed_=0,
00168                         int next_next_=1,
00169                         bool local_random_=true) {
00170     allocate_lattice(ndim_,ndim_,nx_,where_,neighbour_,
00171                      random_seed_, next_next_, local_random_);    
00172   }
00174   void allocate_lattice(int ndim_, 
00175                         int ndir_,
00176                         int nx_[],
00177                         int (*where_)(int*, int, int*)=default_partitioning0,
00178                         void (*neighbour_)(int,int*,int*,int*,int,int*) 
00179                         =torus_topology,
00180                         long random_seed_=0,
00181                         int next_next_=1,
00182                         bool local_random_=true) {
00183     mpi.begin_function("allocate_lattice");
00184     local_random_generator=local_random_;
00185     if(ndim_!=ndir_)
00186       mpi << "It is getting complicated: you have ndim!=ndir\n"; 
00187     deallocate_memory();
00188     // //////////////////////////////////////////////////////////////////
00189     // (*where)(x,nc) must return the processor rank of where x is stored
00190     // (*neghbour)(mu,x_dw,x,x_up,ndim, nx) must fill x_dw and x_up 
00191     //            according with current position x and direction mu
00192     // //////////////////////////////////////////////////////////////////
00193     mpi << "Initializing a mdp_lattice...\n";
00194 
00195     int   mu, nu, rho, process, np;
00196     bool  is_boundary;
00197     long  global_idx, old_idx, new_idx;
00198     ndim=ndim_;
00199     ndir=ndir_;
00200     where=where_;
00201     neighbour=neighbour_;
00202     next_next=next_next_;
00203     nvol=0;
00204     nvol_gl=1;
00205 
00206     mpi << "Lattice dimension: " << nx_[0];
00207     for(mu=1; mu<ndim; mu++) mpi << " x " << nx_[mu];
00208     mpi << '\n';
00209 
00211     // Dynamically allocate some arrays
00213     nx=new int[ndim];
00214 
00215     for(mu=0; mu<ndim; mu++) nvol_gl*=(nx[mu]=nx_[mu]);
00216     int   x[10];
00217     int   x_up[10],x_up_dw[10],x_up_up[10],x_up_up_dw[10],x_up_up_up[10]; 
00218     int   x_dw[10],x_dw_up[10],x_dw_dw[10],x_dw_dw_dw[10],x_dw_dw_up[10];
00219 #if !defined(MDP_NO_LG)
00220     long* local_mdp_sites=new long[nvol_gl];
00221 #else
00222     long lms_tmp=0;
00223     FILE* lms_file=tmpfile();
00224     if(lms_file==0) error("mdp_lattice::generic_lattice()\n"
00225                           "Unable to create temporary lms file");
00226 #endif
00227     // ///////////////////////////////////////////////////////////////////
00228     // Fill table local_mdp_sites with the global coordinate of mdp_sites in ME 
00229     // nvol counts the mdp_sites stored in local_mdp_sites
00230     // ///////////////////////////////////////////////////////////////////
00231     //
00232     // cases of intereset...:
00233     //
00234     // next_next < 0 => only local mdp_sites NEW!
00235     // next_next = 0 => only x+mu and x-mu (wilson) NEW!
00236     // next_next = 1 => only x+-mu, x+-mu+-nu mu!=nu (clover)
00237     // next_next = 2 => only x+-mu, x+-mu+-nu
00238     // next_next = 3 => only x+-mu, x+-mu+-nu, x+-mu+-nu+-rho (asqtad)
00239     // 
00240     // ///////////////////////////////////////////////////////////////////
00241 
00242     for(mu=0; mu<ndim; mu++) x[mu]=0;
00243     do {
00244       global_idx=global_coordinate(x);
00245       if((*where)(x,ndim,nx)==ME) {
00246 #if !defined(MDP_NO_LG)
00247         local_mdp_sites[nvol]=global_idx;
00248 #else
00249         if(fseek(lms_file, nvol*sizeof(long), SEEK_SET)!=0 ||
00250            fwrite(&global_idx, sizeof(long), 1, lms_file)!=1)
00251           error("generic_lattice::allocate_lattice()\n"
00252                 "Unable to write to temporary file");
00253 #endif
00254         nvol++;
00255       } else if(next_next>=0) {
00256         is_boundary=FALSE;
00257         for(mu=0; mu<ndir; mu++) {
00258           (*neighbour)(mu, x_dw, x, x_up, ndim, nx);
00259           if(((*where)(x_up,ndim,nx)>=Nproc) ||
00260              ((*where)(x_dw,ndim,nx)>=Nproc))
00261             error("Incorrect patitioning");
00262           if(((*where)(x_up,ndim,nx)==ME) ||
00263              ((*where)(x_dw,ndim,nx)==ME)) is_boundary=TRUE;
00264           // ////////////////////////////////////////////////////
00265           // cases:
00266           // 1) mu+nu
00267           // 2) mu+nu and mu+mu
00268           // 3) mu+nu+rho and mu+mu+rho and mu+mu+mu
00269           // One may want to optimize case 3. Is was thought for 
00270           // improved staggered fermions. 
00271           // ////////////////////////////////////////////////////
00272           for(nu=0; nu<ndir; nu++)
00273             if((nu!=mu) || (next_next>1)) {
00274               (*neighbour)(nu, x_dw_dw, x_dw, x_dw_up, ndim, nx);      
00275               (*neighbour)(nu, x_up_dw, x_up, x_up_up, ndim, nx);      
00276               if(((*where)(x_up_dw,ndim,nx)>=Nproc) ||
00277                  ((*where)(x_up_up,ndim,nx)>=Nproc) ||
00278                  ((*where)(x_dw_dw,ndim,nx)>=Nproc) ||
00279                  ((*where)(x_dw_up,ndim,nx)>=Nproc)) 
00280                 error("Incorrect patitioning");
00281               if(((*where)(x_up_dw,ndim,nx)==ME) ||
00282                  ((*where)(x_up_up,ndim,nx)==ME) ||
00283                  ((*where)(x_dw_dw,ndim,nx)==ME) ||
00284                  ((*where)(x_dw_up,ndim,nx)==ME)) is_boundary=TRUE;
00285               if(next_next==3)        
00286                 for(rho=0; rho<ndir; rho++) {
00287                   (*neighbour)(rho, x_dw_dw_dw, x_dw_dw, x_dw_dw_up, ndim, nx);      
00288                   (*neighbour)(rho, x_up_up_dw, x_up_up, x_up_up_up, ndim, nx);      
00289                   if(((*where)(x_up_up_up,ndim,nx)>=Nproc) ||
00290                      ((*where)(x_up_up_dw,ndim,nx)>=Nproc) ||
00291                      ((*where)(x_dw_dw_up,ndim,nx)>=Nproc) ||
00292                      ((*where)(x_dw_dw_dw,ndim,nx)>=Nproc)) 
00293                     error("Incorrect patitioning");
00294                   if(((*where)(x_up_up_up,ndim,nx)==ME) ||
00295                      ((*where)(x_up_up_dw,ndim,nx)==ME) ||
00296                      ((*where)(x_dw_dw_up,ndim,nx)==ME) ||
00297                      ((*where)(x_dw_dw_dw,ndim,nx)==ME)) is_boundary=TRUE;
00298                 }
00299             }
00300         }
00301         if(is_boundary==TRUE) {
00302 #if !defined(MDP_NO_LG)
00303           local_mdp_sites[nvol]=global_idx;
00304 #else
00305         if(fseek(lms_file, nvol*sizeof(long), SEEK_SET)!=0 ||
00306            fwrite(&global_idx, sizeof(long), 1, lms_file)!=1)
00307           error("generic_lattice::allocate_lattice()\n"
00308                 "Unable to write to temporary file");
00309 #endif
00310         nvol++;
00311         }
00312       }
00313       x[0]++;
00314       for(mu=0; mu<ndim-1; mu++) if(x[mu]>=nx[mu]) { x[mu]=0; x[mu+1]++; }
00315     } while(x[ndim-1]<nx[ndim-1]);
00316     
00317     // /////////////////////////////////////////////////////////////////
00318     // Dynamically allocate some other arrays
00319     // /////////////////////////////////////////////////////////////////
00320     
00321     dw=new long*[nvol];
00322     up=new long*[nvol];
00323     co=new int*[nvol];
00324     for(new_idx=0; new_idx<nvol; new_idx++) {
00325       dw[new_idx]=new long[ndir];
00326       up[new_idx]=new long[ndir];
00327       co[new_idx]=new int[ndir];
00328     }
00329     gl=new long[nvol];
00330 #if !defined(MDP_NO_LG)
00331     lg=new long[nvol_gl];
00332 #else
00333     lg_file=tmpfile();
00334     if(lg_file==0) error("mdp_lattice::generic_lattice()\n"
00335                          "Unable to create temporary lg file");
00336 #endif
00337     wh=new int[nvol];
00338     for(global_idx=0; global_idx<nvol_gl; global_idx++) 
00339 #if !defined(MDP_NO_LG)
00340       lg[global_idx]=NOWHERE; 
00341 #else
00342     if(fseek(lg_file, global_idx*sizeof(long), SEEK_SET)!=0 ||
00343        fwrite(&NOWHERE, sizeof(long), 1, lg_file)!=1)
00344       error("generic_lattice::allocate_lattice()\n"
00345             "Unable to write to temporary file");
00346 #endif
00347     parity=new int[nvol];
00348     // /////////////////////////////////////////////////////////////////
00349     start[0][0]=stop[0][0]=0;
00350     for(process=0; process<Nproc; process++) {
00351       if(process>0) start[process][0]=stop[process][0]=stop[process-1][1];
00352       for(np=0; np<2; np++) {
00353         if(np>0) start[process][1]=stop[process][1]=stop[process][0];
00354         for(old_idx=0; old_idx<nvol; old_idx++) {
00355 #if !defined(MDP_NO_LG)
00356           global_coordinate(local_mdp_sites[old_idx], x);
00357 #else
00358           if(fseek(lms_file, old_idx*sizeof(long), SEEK_SET)!=0 ||
00359              fread(&lms_tmp, sizeof(long), 1, lms_file)!=1)
00360             error("generic_lattice::allocate_lattice()\n"
00361                   "Unable to read to temporary file");   
00362           global_coordinate(lms_tmp, x);
00363 #endif
00364           if(((*where)(x,ndim,nx)==process) && (compute_parity(x)==np)) {
00365             new_idx=stop[process][np];
00366 #if !defined(MDP_NO_LG)
00367             lg[local_mdp_sites[old_idx]]=new_idx;
00368             gl[new_idx]=local_mdp_sites[old_idx];
00369 #else
00370             
00371             if(fseek(lms_file, old_idx*sizeof(long), SEEK_SET)!=0 ||
00372                fread(&lms_tmp, sizeof(long), 1, lms_file)!=1)
00373             error("generic_lattice::allocate_lattice()\n"
00374                   "Unable to read to temporary file");   
00375             if(fseek(lg_file, lms_tmp*sizeof(long), SEEK_SET)!=0 ||
00376              fwrite(&new_idx, sizeof(long), 1, lg_file)!=1)
00377               error("generic_lattice::allocate_lattice()\n"
00378                     "Unable to write to temporary file");        
00379             gl[new_idx]=lms_tmp;
00380 #endif
00381             wh[new_idx]=process;
00382             parity[new_idx]=compute_parity(x);
00383             stop[process][np]++;
00384           }
00385         }
00386       }
00387     }
00388     // deallcate temporary array
00389 #if !defined(MDP_NO_LG)
00390     delete[] local_mdp_sites;
00391 #else
00392     fclose(lms_file);
00393 #endif
00394     // /////////////////////////
00395     for(new_idx=0; new_idx<nvol; new_idx++) {
00396       global_coordinate(gl[new_idx],x);
00397       for(mu=0; mu<ndim; mu++) co[new_idx][mu]=x[mu];
00398       for(mu=0; mu<ndir; mu++) {
00399         (*neighbour)(mu,x_dw,x,x_up,ndim,nx);
00400         if(wh[new_idx]==ME) {
00401           dw[new_idx][mu]=local(global_coordinate(x_dw));
00402           up[new_idx][mu]=local(global_coordinate(x_up));
00403         } else {
00404           if(local(global_coordinate(x_dw))!=NOWHERE) 
00405             dw[new_idx][mu]=local(global_coordinate(x_dw));
00406           else dw[new_idx][mu]=NOWHERE;
00407           if(local(global_coordinate(x_up))!=NOWHERE)
00408             up[new_idx][mu]=local(global_coordinate(x_up));
00409           else up[new_idx][mu]=NOWHERE;
00410         }
00411       }
00412     }
00413     nvol_in=stop[ME][1]-start[ME][0];
00414     mpi << "Communicating...\n";
00415     communicate_results_to_all_processes();
00416     mpi << "Initializing random per mdp_site...\n";
00417     if(mdp_random_seed_filename && random_seed_==0) {
00418       if(ME==0) {
00419         FILE *fp=fopen(mdp_random_seed_filename, "r");
00420         if(fp!=0) {
00421           if(fread(&random_seed_, sizeof(random_seed_),1,fp)!=1) random_seed_=0;
00422           fclose(fp);
00423         };
00424         fp=fopen(mdp_random_seed_filename, "w");
00425         if(fp!=0) {
00426           random_seed_+=1;
00427           fwrite(&random_seed_, sizeof(random_seed_),1,fp);
00428           random_seed_-=1;
00429           fclose(fp);
00430         }
00431         mpi << "Reading from file " << mdp_random_seed_filename << " lattice().random_seed=" << random_seed_ << '\n';
00432         mpi << "Writing to   file " << mdp_random_seed_filename << " lattice().random_seed=" << random_seed_+1 << '\n';
00433       }
00434       mpi.broadcast(random_seed_,0);
00435     } else {
00436       mpi << "Adopting random_seed=" << random_seed_ << '\n';
00437     }
00438     initialize_random(random_seed_);
00439     mpi << "Lattice created.\n";
00440     mpi.end_function("allocate_lattice");
00441   }
00442   // ////////////////////////////////////////////////////
00443   // deallocate all the dynamically alocated arrays
00444   // ////////////////////////////////////////////////////
00445   virtual ~mdp_lattice() {
00446     deallocate_memory();
00447   }
00449   void deallocate_memory() {
00450     if(nvol==0) return;
00451     int process, new_idx;
00452     delete[] nx;
00453     for(process=0; process<Nproc; process++) if(process!=ME) {
00454       if(len_to_send[process]!=0) delete[] to_send[process];
00455     }
00456     for(new_idx=0; new_idx<nvol; new_idx++) {
00457       delete[] dw[new_idx];
00458       delete[] up[new_idx];
00459       delete[] co[new_idx];
00460     }
00461     delete[] dw;
00462     delete[] up;
00463     delete[] co;
00464     delete[] wh;
00465     delete[] gl;
00466 #if !defined(MDP_NO_LG)
00467     delete[] lg;
00468 #else
00469     fclose(lg_file);
00470 #endif
00471     delete[] parity;
00472     delete[] random_obj;
00473   }
00474   // ////////////////////////////////////////////////////
00475   // initialize random number generator for each local mdp_site
00476   // ////////////////////////////////////////////////////
00477   void initialize_random(long random_seed_=0) {
00478     random_seed=random_seed_;
00479     if(local_random_generator) {
00480       mdp << "Using a local random generator\n";
00481       if(random_obj!=0) delete[] random_obj;
00482       random_obj=new mdp_prng[nvol_in];
00483       for(long idx=0; idx<nvol_in; idx++)
00484         random_obj[idx].initialize(gl[idx+start[ME][0]]+random_seed);
00485     }
00486   }
00487   inline mdp_prng &random(mdp_site);
00488   // //////////////////////////////////
00489   // functions for external access ... 
00490   // to be used to access member variables
00491   // /////////////////////////////////
00492 
00494   inline int n_dimensions() const {
00495     return ndim;
00496   }
00498   inline int n_directions() const {
00499     return ndir;
00500   }
00502   inline long size() const {
00503     return nvol_gl;
00504   }
00506   inline long size(const int mu) const {
00507     return nx[mu];
00508   }
00510   inline long local_volume() const {
00511     return nvol_in;
00512   }
00514   inline long global_volume() const {
00515     return nvol_gl;
00516   }
00517   inline long move_up(const long idx, const int mu) const {
00518     return up[idx][mu];
00519   }
00520   inline long move_down(const long idx, const int mu) const {
00521     return up[idx][mu];
00522 }
00523   inline long local(long idx) const {
00524 #if !defined(MDP_NO_LG)
00525     return lg[idx];
00526 #else
00527     long lg_tmp;
00528     if(fseek(lg_file, idx*sizeof(long), SEEK_SET)!=0 ||
00529        fread(&lg_tmp, sizeof(long), 1, lg_file)!=1)
00530       error("generic_lattice::allocate_lattice()\n"
00531             "Unable to read to temporary file");         
00532     return lg_tmp;
00533 #endif
00534   }
00535   inline long global(long idx) const {
00536     return gl[idx];
00537   }
00538   inline int site_parity(const long idx) const {
00539     return parity[idx];
00540   }
00541   inline long start_index(const int process, int p=EVENODD) const {
00542     if(p==EVENODD) p=0;
00543     return start[process][p];
00544   }
00545   inline long stop_index(const int process, int p=EVENODD) const {
00546     if(p==EVENODD) p=1;
00547     return stop[process][p];
00548   }
00549 };

Generated on Sun Feb 27 15:12:21 2005 by  doxygen 1.4.1