fermiqcd: /Users/mdipierro/fermiqcd/development/Libraries/fermiqcd_gauge

00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 /*
00014 void mul_OUU(mdp_complex &A, mdp_complex &B, mdp_complex &C,
00015              int nc=3, bool init=false) {
00016   nc=B.rows();
00017   for(int i=0; i<nc; i++)
00018     for(int j=0; j<nc; j++) {
00019       if(init)(&A)[i*nc+j]=0.0;
00020       for(int k=0; k<nc; k++)
00021         (&A)[i*nc+j]+=(&B)[i*nc+k]*(&C)[k*nc+j];      
00022     }
00023 }
00024 */
00025 
00027 class gauge_stats {};
00028 
00029 
00042 class WilsonGaugeAction {
00043  public:
00044   static void heatbath_SU2(mdp_prng &random,
00045                            mdp_real beta_eff, mdp_complex *a) {
00046     mdp_real e0,e1,e2,e3, dk, p0;
00047     mdp_real r1,r2,r3,r4;
00048     mdp_real a0,a1,a2,a3;
00049     mdp_complex u0,u1,u2,u3;
00050     mdp_complex v0,v1,v2,v3;
00051     mdp_real delta, phi, sin_alpha, sin_theta, cos_theta;
00052     e0=real(a[0])+real(a[3]);
00053     e1=imag(a[1])+imag(a[2]);
00054     e2=real(a[1])-real(a[2]);
00055     e3=imag(a[0])-imag(a[3]);
00056     dk=sqrt(e0*e0+e1*e1+e2*e2+e3*e3);
00057     p0=(dk*beta_eff);
00058     u0=mdp_complex(e0/dk,-e3/dk);
00059     u2=mdp_complex(e2/dk,-e1/dk);
00060     u1=mdp_complex(-e2/dk,-e1/dk);
00061     u3=mdp_complex(e0/dk,e3/dk);
00062     
00063     if(beta_eff<=0.0) 
00064       error("fermiqcd_gauge_algorithms/heatbath_SU2: beta is zero");
00065     
00066     do {
00067       do; while ((r1 = random.plain()) < 0.0001);
00068       r1 = -log(r1)/p0;
00069       do; while ((r2 = random.plain()) < 0.0001);
00070       r2 = -log(r2)/p0;
00071       r3 = cos(2.0*Pi*random.plain());
00072       r3 = r3*r3;
00073       delta = r2+r1*r3;
00074       r4=random.plain();
00075     } while (r4*r4 > (1.0-0.5*delta));
00076     a0=1.0-delta;
00077     cos_theta=2.0*random.plain()-1.0;
00078     sin_theta=sqrt(1.0-cos_theta*cos_theta);
00079     sin_alpha=sqrt(1-a0*a0);
00080     phi=2.0*Pi*random.plain();
00081     a1=sin_alpha*sin_theta*cos(phi);
00082     a2=sin_alpha*sin_theta*sin(phi);
00083     a3=sin_alpha*cos_theta;
00084     v0=mdp_complex(a0,a3);
00085     v1=mdp_complex(a2,a1);
00086     v2=mdp_complex(-a2,a1);
00087     v3=mdp_complex(a0,-a3);
00088     a[0]=v0*u0+v1*u2;
00089     a[1]=v0*u1+v1*u3;
00090     a[2]=v2*u0+v3*u2;
00091     a[3]=v2*u1+v3*u3;
00092   }
00093 
00094   /*
00095   static void contract_staples(mdp_matrix &M, gauge_field &U, mdp_site &U, int mu){
00096     int nc=U.nc;
00097     mdp_matrix B(nc,nc);
00098     mdp_site y=x+mu;
00099     mdp_site z;
00100     for(int nu=0; nu<U.ndim; nu++) if(nu!=mu) {
00101         z=y+nu;
00102         mul_OUH(M(0,0),U(y,nu,0,0),U(x+nu,mu,0,0));
00103         mul_OUH(M(0,0),U(y,nu,0,0),U(x+nu,mu,0,0));
00104       }
00105     
00106   }
00107   */
00108   static gauge_stats heatbath(gauge_field &U, 
00109                               coefficients &coeff, 
00110                               int n_iter=1) {
00111     begin_function("WilsonGaugeAction__heatbath");
00112     if(U.nc==1) error("fermiqcd_gauge_algorithms/heatbath(): U(1)? (use metropolis)");
00113     gauge_stats stats;
00114     mdp_real beta, zeta;
00115     if(coeff.has_key("beta")) beta=coeff["beta"]; 
00116     else error("beta undeclared");
00117     if(coeff.has_key("zeta")) zeta=coeff["zeta"]; 
00118     else zeta=1;
00119 
00120     int i,j,k,iter,mu,parity;
00121     mdp_matrix M(U.nc,U.nc);
00122     mdp_complex a[4], tmpUik;
00123     site x(U.lattice());
00124     double time=mpi.time();
00125 
00126     mdp << coeff;
00127 
00128     for(iter=0; iter<n_iter; iter++) 
00129       for(parity=0; parity<2; parity++) 
00130         for(mu=0; mu<U.ndim; mu++) {
00131           forallsitesofparity(x,parity) {
00132             for(i=0; i<U.nc-1; i++)
00133               for(j=i+1; j<U.nc; j++) {
00134                 if(zeta==1)    M=U(x,mu)*staple_H(U,x,mu);
00135                 else if(mu==0) M=zeta*U(x,0)*staple_H(U,x,0);
00136                 else           M=((mdp_real) 1.0/zeta)*U(x,mu)*staple_H(U,x,mu);
00137                 a[0]=M(i,i); 
00138                 a[1]=M(i,j);
00139                 a[2]=M(j,i);
00140                 a[3]=M(j,j);
00141                 heatbath_SU2(U.lattice().random(x),beta/U.nc,a);
00142                 for(k=0; k<U.nc; k++) {
00143                   tmpUik=a[0]*U(x,mu,i,k)+a[1]*U(x,mu,j,k);
00144                   U(x,mu,j,k)=a[2]*U(x,mu,i,k)+a[3]*U(x,mu,j,k);
00145                   U(x,mu,i,k)=tmpUik;
00146                 }
00147               }
00148           }
00149           // The next command does all the communications!
00150           U.update(parity, mu, U.nc*U.nc);
00151         }
00152     mdp << "\t<stats>\n\t\t<time>" << mpi.time()-time << "</time>\n\t</stats>\n";
00153     end_function("WilsonGaugeAction__heatbath");
00154     return stats;
00155   }
00156 };
00157 
00186 class ImprovedGaugeAction : public WilsonGaugeAction {
00187  private:
00188   static mdp_matrix rectangles_0i_H(gauge_field &U, site x, int mu) {
00189     mdp_matrix tmp(U.nc,U.nc);
00190     mdp_matrix b1(U.nc, U.nc);
00191     mdp_matrix b2(U.nc, U.nc);
00192     mdp_matrix b3(U.nc, U.nc);
00193     site y0(U.lattice());
00194     site y1(U.lattice());
00195     site y2(U.lattice());
00196     int nu;
00197     tmp=0;
00198     if(mu==0) {
00199       for(nu=1; nu<U.ndim; nu++) {
00200         y0=x+mu;
00201         y1=y0+nu;
00202         tmp+=U(y0,nu)*U(y1,nu)*hermitian(U(x,nu)*U(x+nu,nu)*U((x+nu)+nu,mu));
00203         y0=(x-nu)-nu;
00204         y1=y0+mu;
00205         y2=y1+nu;
00206         tmp+=hermitian(U(y0,mu)*U(y1,nu)*U(y2,nu))*U(y0,nu)*U(x-nu,nu);
00207       }
00208     } else {
00209       nu=0;
00210       y0=(x-mu)+nu;
00211       tmp+=U(x+mu,nu)*hermitian(U(x-mu,nu)*U(y0,mu)*U(x+nu,mu))*U(x-mu,mu);
00212       y0=x-mu;
00213       y1=y0-nu;
00214       y2=y1+mu;
00215       tmp+=hermitian(U(y1,mu)*U(y2,mu)*U(y2+mu,nu))*U(y1,nu)*U(y0,mu);
00216       y0=x+mu;
00217       y1=y0+nu;
00218       y2=y1-mu;
00219       tmp+=U(y0,mu)*U(y0+mu,nu)*hermitian(U(x,nu)*U(y2,mu)*U(y1,mu));
00220       y0=((x+mu)+mu)-nu;
00221       y1=y0-mu;
00222       y2=y1-mu;
00223       tmp+=U(x+mu,mu)*hermitian(U(y2,mu)*U(y1,mu)*U(y0,nu))*U(y2,nu);
00224       
00225     }
00226     prepare(tmp);
00227     return tmp;
00228   }
00229   
00230 // if min_nu==0 then rectangles_ij computes all 6 rectanges
00231   
00232   static mdp_matrix rectangles_ij_H(gauge_field &U, site x, int mu, int min_nu=1) {
00233     mdp_matrix tmp(U.nc,U.nc);
00234     mdp_matrix b1(U.nc, U.nc);
00235     mdp_matrix b2(U.nc, U.nc);
00236     mdp_matrix b3(U.nc, U.nc);
00237     site y0(U.lattice());
00238     site y1(U.lattice());
00239     site y2(U.lattice());
00240     int nu;
00241     for(nu=min_nu; nu<U.ndim; nu++) if(nu!=mu) {    
00242       y0=x+mu;
00243       y1=y0+nu;
00244       tmp+=U(y0,nu)*U(y1,nu)*hermitian(U(x,nu)*U(x+nu,nu)*U((x+nu)+nu,mu));
00245       
00246       y0=(x-nu)-nu;
00247       y1=y0+mu;
00248       y2=y1+nu;
00249       tmp+=hermitian(U(y0,mu)*U(y1,nu)*U(y2,nu))*U(y0,nu)*U(x-nu,nu);
00250       
00251       y0=(x-mu)+nu;
00252       tmp+=U(x+mu,nu)*hermitian(U(x-mu,nu)*U(y0,mu)*U(x+nu,mu))*U(x-mu,mu);
00253       
00254       y0=x-mu;
00255       y1=y0-nu;
00256       y2=y1+mu;
00257       tmp+=hermitian(U(y1,mu)*U(y2,mu)*U(y2+mu,nu))*U(y1,nu)*U(y0,mu);
00258 
00259       y0=x+mu;
00260       y1=y0+nu;
00261       y2=y1-mu;
00262       tmp+=U(y0,mu)*U(y0+mu,nu)*hermitian(U(x,nu)*U(y2,mu)*U(y1,mu));
00263 
00264       y0=((x+mu)+mu)-nu;
00265       y1=y0-mu;
00266       y2=y1-mu;
00267       tmp+=U(x+mu,mu)*hermitian(U(y2,mu)*U(y1,mu)*U(y0,nu))*U(y2,nu);
00268     }
00269     prepare(tmp);
00270     return tmp;
00271   }
00272 
00273   // //////////////////////////////////////////////////////
00274   // this is slow but should make the chair correcly ...
00275   // see: hep-lat/0712010
00276   // //////////////////////////////////////////////////////
00277 
00278   static mdp_matrix chair_H(gauge_field &U, site x, int mu) {
00279     int ndim=U.ndim;
00280     int nu, rho;
00281     mdp_matrix tmp(U.nc, U.nc);
00282     mdp_matrix b1(U.nc, U.nc);
00283     mdp_matrix b2(U.nc, U.nc);
00284     mdp_matrix b3(U.nc, U.nc);
00285     site y1(U.lattice());
00286     site y2(U.lattice());
00287     site y3(U.lattice());
00288     site y4(U.lattice());
00289     site y5(U.lattice());
00290     tmp=0;
00291     for(nu=0; nu<ndim; nu++) if(nu!=mu)
00292       for(rho=0; rho<ndim; rho++) if((rho!=nu) && (rho!=mu)) { 
00293         y1=x+mu;
00294         y2=y1+nu;
00295         y3=y2+rho;
00296         y4=y3-mu;
00297         y5=y4-nu;
00298         tmp+=U(y1,nu)*U(y2,rho)*hermitian(U(x,rho)*U(y5,nu)*U(y4,mu));
00299         y1=x+mu;
00300         y2=y1-nu;
00301         y3=y2+rho;
00302         y4=y3-mu;
00303         y5=y4+nu;
00304         tmp+=hermitian(U(y2,nu))*U(y2,rho)
00305           *hermitian(U(y4,mu))*U(y4,nu)*hermitian(U(x,rho));
00306         
00307         y1=x+mu;
00308         y2=y1+nu;
00309         y3=y2-rho;
00310         y4=y3-mu;
00311         y5=y4-nu;
00312         tmp+=U(y1,nu)*hermitian(U(y5,nu)*U(y4,mu)*U(y3,rho))*U(y5,rho);
00313         
00314         y1=x+mu;
00315         y2=y1-nu;
00316         y3=y2-rho;
00317         y4=y3-mu;
00318         y5=y4+nu;
00319         tmp+=hermitian(U(y4,mu)*U(y3,rho)*U(y2,nu))*U(y4,nu)*U(y5,rho);
00320 
00321       }
00322     return(tmp);
00323   }
00324   
00325   // ////////////////////////////////////////////////////////////////////
00326   // new_heatbath uses an improved gauge action!
00327   // both isotropic (param.zeta=1) and anisotropic
00328   // ////////////////////////////////////////////////////////////////////
00329 
00330   enum { iGauge_min=4 };
00331 
00332   static int strange_mapping(site &x) {
00333     int mu, type=0;
00334     for(mu=0; mu<x.lattice().ndim; mu++) type+=(int) pow((float) iGauge_min,mu)*(x(mu) % iGauge_min);
00335     return type;
00336   }
00337   
00338  public:
00339   static gauge_stats heatbath(gauge_field &U,
00340                               coefficients &coeff,
00341                               int n_iter=1,
00342                               string model="MILC") { 
00343     
00344     begin_function("ImprovedGaugeAction__heatbath");
00345 
00346     gauge_stats stats;
00347     mdp_real beta, zeta, u_s, u_t;
00348 
00349     if(coeff.has_key("beta")) beta=coeff["beta"]; else error("beta undeclared");
00350     if(coeff.has_key("zeta")) zeta=coeff["zeta"]; else zeta=1;
00351     if(coeff.has_key("u_t"))  u_t=coeff["u_t"];   else u_t=1;
00352     if(coeff.has_key("u_s"))  u_s=coeff["u_s"];   else u_s=1;
00353     
00354     // if(Nproc!=1)  error("improved_heatbath() does not work in parallel!");
00355 
00356     if(U.ndim!=4) error("fermiqcd_gauge_algorithms/improved_heatbath(): ndim!=4 (use heatbath instead)");
00357 
00358     /*
00359       if((U.lattice().size(0) % 3!=0) || (U.lattice().size(1) % 3!=0) ||
00360       (U.lattice().size(2) % 3!=0) || (U.lattice().size(3) % 3!=0))
00361       error("lattice is not divisible by 3");
00362     */
00363 
00364     if(U.nc==1)   error("fermiqcd_gauge_algorithms/improved_heatbath(): U(1)? (use metropolis instead)");
00365     int nc=U.nc;
00366     int ndim=U.ndim;
00367     int i,j,k,iter,mu,type;
00368     mdp_matrix M;
00369     site x(U.lattice());
00370     double time=mpi.time();
00371     mdp_complex a[4],tmpUik;
00372     mdp_real alpha_s;
00373     mdp_real c_tp=0, c_tr=0, c_sp=0, c_sr=0, c_p=0, c_r=0, c_c=0;
00374     
00375     if(model=="Morningstar") {
00376       
00377       c_tp=4.0*zeta/(3.0*pow((double)u_s*u_t,(double)2.0));
00378       c_tr=-1.0*zeta/(12.0*pow((double)u_s,(double)4.0)*pow((double)u_t,(double)2.0));
00379       c_sp=5.0/(3.0*zeta*pow((double)u_s,(double)4.0));
00380       c_sr=-1.0/(12.0*zeta*pow((double)u_s,(double)6.0));
00381       
00382       c_p=1.0*pow((double)u_s,(double)-4.0);
00383       c_r=-0.05*pow((double)u_s,(double)-6.0);
00384       c_c=0;
00385       
00386     } else if(model=="MILC") {
00387       
00388       if(zeta!=1) 
00389         error("fermiqcd_gauge_algorithms/improved_heatbath: zeta!=1");
00390       
00391       alpha_s=-4.0*log(u_s)/3.0684;
00392       c_p=1.0;
00393       c_r=-0.05*pow((double)u_s,(double)-2.0)*(1.0+0.4805*alpha_s);;
00394       c_c=-1.00*pow((double)u_s,(double)-2.0)*(0.03325*alpha_s);;
00395       
00396     } else {
00397       mdp << "Using default non-improved action" << '\n';
00398       stats=WilsonGaugeAction::heatbath(U,coeff,n_iter);
00399       end_function("ImprovedGaugeAction__heatbath");
00400       return stats;
00401     }
00402     
00403     mdp << coeff;
00404     
00405     for(iter=0; iter<n_iter; iter++)
00406       for(type=0; type<(int) pow((float) iGauge_min, U.ndim); type++) {
00407         forallsites(x) {
00408           if(strange_mapping(x)==type) {
00409             for(mu=0; mu<ndim; mu++) 
00410               for(i=0; i<nc-1; i++)
00411                 for(j=i+1; j<nc; j++) { 
00412                   if(zeta!=1) {
00413                     // anisotropic San Diego action
00414                     if(mu==0) M=U(x,mu)*(c_tp*staple_0i_H(U,x,0)+
00415                                          c_tr*rectangles_0i_H(U,x,0));
00416                     else      M=U(x,mu)*(c_sp*staple_ij_H(U,x,mu)+
00417                                          c_tp*staple_0i_H(U,x,mu)+
00418                                          c_sr*rectangles_ij_H(U,x,mu)+
00419                                          c_tr*rectangles_0i_H(U,x,mu));
00420                   } else if(c_c==0) {
00421                     // isotropic San Diego action
00422                     M=U(x,mu)*(c_p*staple_H(U,x,mu)+
00423                                c_r*rectangles_ij_H(U,x,mu,0));
00424                   } else {
00425                     // fully O(a^2) improved isotropic MILC action
00426                     M=U(x,mu)*(c_p*staple_H(U,x,mu)+
00427                                c_r*rectangles_ij_H(U,x,mu,0)+
00428                                c_c*chair_H(U,x,mu));
00429                   }
00430                   a[0]=M(i,i); 
00431                   a[1]=M(i,j);
00432                   a[2]=M(j,i);
00433                   a[3]=M(j,j);
00434                   heatbath_SU2(U.lattice().random(x),beta/U.nc,a);
00435                   for(k=0; k<U.nc; k++) {
00436                     tmpUik=a[0]*U(x,mu,i,k)+a[1]*U(x,mu,j,k);
00437                     U(x,mu,j,k)=a[2]*U(x,mu,i,k)+a[3]*U(x,mu,j,k);
00438                     U(x,mu,i,k)=tmpUik;
00439                   }
00440                 }
00441           }
00442         }
00443         U.update();
00444       }
00445     mdp << "\t<stats>\n\t\t<time>" << mpi.time()-time << "</time>\n\t</stats>\n";
00446     end_function("ImprovedGaugeAction__heatbath");
00447     return stats;
00448   }
00449 };
00450
/Users/mdipierro/fermiqcd/development/Libraries/fermiqcd_gauge_actions.h
