fermiqcd_gauge_actions_sse2.h Source File

00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00042 class ImprovedGaugeActionSSE2 : public WilsonGaugeAction {
00043  private:
00044 
00045 #if !defined(SSE2) && !defined(USE_DOUBLE_PRECISION)
00046 
00047   static mdp_matrix rectangles_0i_H(gauge_field &U, site x, int mu) {
00048     mdp_matrix tmp(3,3);
00049     mdp_matrix b1(3, 3);
00050     mdp_matrix b2(3, 3);
00051     mdp_matrix b3(3, 3);
00052     site y0(U.lattice());
00053     site y1(U.lattice());
00054     site y2(U.lattice());
00055     int nu;
00056     tmp=0;
00057     if(mu==0) {
00058       for(nu=1; nu<U.ndim; nu++) {
00059         y0=x+mu;
00060         y1=y0+nu;
00061         tmp+=U(y0,nu)*U(y1,nu)*hermitian(U(x,nu)*U(x+nu,nu)*U((x+nu)+nu,mu));
00062         y0=(x-nu)-nu;
00063         y1=y0+mu;
00064         y2=y1+nu;
00065         tmp+=hermitian(U(y0,mu)*U(y1,nu)*U(y2,nu))*U(y0,nu)*U(x-nu,nu);
00066       }
00067     } else {
00068       nu=0;
00069       y0=(x-mu)+nu;
00070       tmp+=U(x+mu,nu)*hermitian(U(x-mu,nu)*U(y0,mu)*U(x+nu,mu))*U(x-mu,mu);
00071       y0=x-mu;
00072       y1=y0-nu;
00073       y2=y1+mu;
00074       tmp+=hermitian(U(y1,mu)*U(y2,mu)*U(y2+mu,nu))*U(y1,nu)*U(y0,mu);
00075       y0=x+mu;
00076       y1=y0+nu;
00077       y2=y1-mu;
00078       tmp+=U(y0,mu)*U(y0+mu,nu)*hermitian(U(x,nu)*U(y2,mu)*U(y1,mu));
00079       y0=((x+mu)+mu)-nu;
00080       y1=y0-mu;
00081       y2=y1-mu;
00082       tmp+=U(x+mu,mu)*hermitian(U(y2,mu)*U(y1,mu)*U(y0,nu))*U(y2,nu);
00083 
00084     }
00085 
00086     return tmp;
00087   }
00088 
00089   // if min_nu==0 then rectangles_ij computes all 6 rectanges
00090 
00091   static mdp_matrix rectangles_ij_H(gauge_field &U, site x, int mu, int min_nu=1) {
00092     mdp_matrix tmp(3,3);
00093     mdp_matrix b1(3, 3);
00094     mdp_matrix b2(3, 3);
00095     mdp_matrix b3(3, 3);
00096     site y0(U.lattice());
00097     site y1(U.lattice());
00098     site y2(U.lattice());
00099     int nu;
00100     for(nu=min_nu; nu<U.ndim; nu++) if(nu!=mu) {    
00101       y0=x+mu;
00102       y1=y0+nu;
00103       tmp+=U(y0,nu)*U(y1,nu)*hermitian(U(x,nu)*U(x+nu,nu)*U((x+nu)+nu,mu));
00104       
00105       y0=(x-nu)-nu;
00106       y1=y0+mu;
00107       y2=y1+nu;
00108       tmp+=hermitian(U(y0,mu)*U(y1,nu)*U(y2,nu))*U(y0,nu)*U(x-nu,nu);
00109       
00110       y0=(x-mu)+nu;
00111       tmp+=U(x+mu,nu)*hermitian(U(x-mu,nu)*U(y0,mu)*U(x+nu,mu))*U(x-mu,mu);
00112       
00113       y0=x-mu;
00114       y1=y0-nu;
00115       y2=y1+mu;
00116       tmp+=hermitian(U(y1,mu)*U(y2,mu)*U(y2+mu,nu))*U(y1,nu)*U(y0,mu);
00117 
00118       y0=x+mu;
00119       y1=y0+nu;
00120       y2=y1-mu;
00121       tmp+=U(y0,mu)*U(y0+mu,nu)*hermitian(U(x,nu)*U(y2,mu)*U(y1,mu));
00122 
00123       y0=((x+mu)+mu)-nu;
00124       y1=y0-mu;
00125       y2=y1-mu;
00126       tmp+=U(x+mu,mu)*hermitian(U(y2,mu)*U(y1,mu)*U(y0,nu))*U(y2,nu);
00127     }
00128     return tmp;
00129   }
00130 
00131   // //////////////////////////////////////////////////////
00132   // this is slow but should make the chair correcly ...
00133   // see: hep-lat/0712010
00134   // //////////////////////////////////////////////////////
00135 
00136   static mdp_matrix chair_H(gauge_field &U, site x, int mu) {
00137     int ndim=U.ndim;
00138     int nu, rho;
00139     mdp_matrix tmp(3, 3);
00140     mdp_matrix b1(3, 3);
00141     mdp_matrix b2(3, 3);
00142     mdp_matrix b3(3, 3);
00143     site y1(U.lattice());
00144     site y2(U.lattice());
00145     site y3(U.lattice());
00146     site y4(U.lattice());
00147     site y5(U.lattice());
00148     tmp=0;
00149     for(nu=0; nu<ndim; nu++) if(nu!=mu)
00150       for(rho=0; rho<ndim; rho++) if((rho!=nu) && (rho!=mu)) { 
00151         y1=x+mu;
00152         y2=y1+nu;
00153         y3=y2+rho;
00154         y4=y3-mu;
00155         y5=y4-nu;
00156         tmp+=U(y1,nu)*U(y2,rho)*hermitian(U(x,rho)*U(y5,nu)*U(y4,mu));
00157         y1=x+mu;
00158         y2=y1-nu;
00159         y3=y2+rho;
00160         y4=y3-mu;
00161         y5=y4+nu;
00162         tmp+=hermitian(U(y2,nu))*U(y2,rho)
00163           *hermitian(U(y4,mu))*U(y4,nu)*hermitian(U(x,rho));
00164         
00165         y1=x+mu;
00166         y2=y1+nu;
00167         y3=y2-rho;
00168         y4=y3-mu;
00169         y5=y4-nu;
00170         tmp+=U(y1,nu)*hermitian(U(y5,nu)*U(y4,mu)*U(y3,rho))*U(y5,rho);
00171         
00172         y1=x+mu;
00173         y2=y1-nu;
00174         y3=y2-rho;
00175         y4=y3-mu;
00176         y5=y4+nu;
00177         tmp+=hermitian(U(y4,mu)*U(y3,rho)*U(y2,nu))*U(y4,nu)*U(y5,rho);
00178         
00179       }
00180     return(tmp);
00181   }
00182   
00183 #else
00184   
00185   static mdp_matrix rectangles_0i_H(gauge_field &U, site x, int mu) {
00186     int nc=3;
00187     mdp_matrix tmp(nc,nc);
00188     mdp_matrix b1(nc, nc);
00189     mdp_matrix b2(nc, nc);
00190     mdp_matrix b3(nc, nc);
00191     site y0(U.lattice());
00192     site y1(U.lattice());
00193     site y2(U.lattice());
00194     int nu;
00195     tmp=0;
00196     if(mu==0) {
00197       for(nu=1; nu<U.ndim; nu++) {
00198         y0=x+mu;
00199         y1=y0+nu;
00200         _sse_mulABC_set_333(&U(y0,nu,0,0),&U(y1,nu,0,0),&b1(0,0));
00201         _sse_mulABHC_set_333(&b1(0,0),&U((x+nu)+nu,mu,0,0),&b2(0,0));
00202         _sse_mulABHC_set_333(&b2(0,0),&U(x+nu,nu,0,0),&b3(0,0));
00203         _sse_mulABHC_add_333(&b3(0,0),&U(x,nu,0,0),&tmp(0,0));
00204         y0=(x-nu)-nu;
00205         y1=y0+mu;
00206         y2=y1+nu;
00207         _sse_mulABC_set_333(&U(y0,nu,0,0),&U(x-nu,nu,0,0),&b1(0,0));
00208         _sse_mulAHBC_set_333(&U(y0,mu,0,0),&b1(0,0),&b2(0,0));
00209         _sse_mulAHBC_set_333(&U(y1,nu,0,0),&b2(0,0),&b3(0,0));
00210         _sse_mulAHBC_add_333(&U(y2,nu,0,0),&b3(0,0),&tmp(0,0));
00211       }
00212     } else {
00213       nu=0;
00214       y0=(x-mu)+nu;
00215       _sse_mulABHC_set_333(&U(x+mu,nu,0,0),&U(x+nu,mu,0,0),&b1(0,0));
00216       _sse_mulABHC_set_333(&b1(0,0),&U(y0,mu,0,0),&b2(0,0));
00217       _sse_mulABHC_set_333(&b2(0,0),&U(x-mu,nu,0,0),&b3(0,0));
00218       _sse_mulABC_add_333(&b3(0,0),&U(x-mu,mu,0,0),&tmp(0,0));
00219       y0=x-mu;
00220       y1=y0-nu;
00221       y2=y1+mu;
00222       _sse_mulABC_set_333(&U(y1,nu,0,0),&U(y0,mu,0,0),&b1(0,0));
00223       _sse_mulAHBC_set_333(&U(y1,mu,0,0),&b1(0,0),&b2(0,0));
00224       _sse_mulAHBC_set_333(&U(y2,mu,0,0),&b2(0,0),&b3(0,0));
00225       _sse_mulAHBC_add_333(&U(y2+mu,nu,0,0),&b3(0,0),&tmp(0,0));
00226       y0=x+mu;
00227       y1=y0+nu;
00228       y2=y1-mu;
00229       _sse_mulABC_set_333(&U(y0,mu,0,0),&U(y0+mu,nu,0,0),&b1(0,0));
00230       _sse_mulABHC_set_333(&b1(0,0),&U(y1,mu,0,0),&b2(0,0));
00231       _sse_mulABHC_set_333(&b2(0,0),&U(y2,mu,0,0),&b3(0,0));
00232       _sse_mulABHC_add_333(&b3(0,0),&U(x,nu,0,0),&tmp(0,0));
00233       y0=((x+mu)+mu)-nu;
00234       y1=y0-mu;
00235       y2=y1-mu;
00236       _sse_mulABHC_set_333(&U(x+mu,mu,0,0),&U(y0,nu,0,0),&b1(0,0));
00237       _sse_mulABHC_set_333(&b1(0,0),&U(y1,mu,0,0),&b2(0,0));
00238       _sse_mulABHC_set_333(&b2(0,0),&U(y2,mu,0,0),&b3(0,0));
00239       _sse_mulABC_add_333(&b3(0,0),&U(y2,nu,0,0),&tmp(0,0));
00240     }
00241     return tmp;
00242   }
00243   
00244   // if min_nu==0 then rectangles_ij computes all 6 rectanges
00245   
00246   static mdp_matrix rectangles_ij_H(gauge_field &U, site x, int mu, int min_nu=1) {
00247     int nc=3;
00248     mdp_matrix tmp(nc,nc);
00249     mdp_matrix b1(nc, nc);
00250     mdp_matrix b2(nc, nc);
00251     mdp_matrix b3(nc, nc);
00252     site y0(U.lattice());
00253     site y1(U.lattice());
00254     site y2(U.lattice());
00255     int nu;
00256     for(nu=min_nu; nu<U.ndim; nu++) if(nu!=mu) {    
00257       y0=x+mu;
00258       y1=y0+nu;
00259       _sse_mulABC_set_333(&U(y0,nu,0,0),&U(y1,nu,0,0),&b1(0,0));
00260       _sse_mulABHC_set_333(&b1(0,0),&U((x+nu)+nu,mu,0,0),&b2(0,0));
00261       _sse_mulABHC_set_333(&b2(0,0),&U(x+nu,nu,0,0),&b3(0,0));
00262       _sse_mulABHC_add_333(&b3(0,0),&U(x,nu,0,0),&tmp(0,0));
00263       y0=(x-nu)-nu;
00264       y1=y0+mu;
00265       y2=y1+nu;
00266       _sse_mulABC_set_333(&U(y0,nu,0,0),&U(x-nu,nu,0,0),&b1(0,0));
00267       _sse_mulAHBC_set_333(&U(y0,mu,0,0),&b1(0,0),&b2(0,0));
00268       _sse_mulAHBC_set_333(&U(y1,nu,0,0),&b2(0,0),&b3(0,0));
00269       _sse_mulAHBC_add_333(&U(y2,nu,0,0),&b3(0,0),&tmp(0,0));
00270       y0=(x-mu)+nu;
00271       _sse_mulABHC_set_333(&U(x+mu,nu,0,0),&U(x+nu,mu,0,0),&b1(0,0));
00272       _sse_mulABHC_set_333(&b1(0,0),&U(y0,mu,0,0),&b2(0,0));
00273       _sse_mulABHC_set_333(&b2(0,0),&U(x-mu,nu,0,0),&b3(0,0));
00274       _sse_mulABC_add_333(&b3(0,0),&U(x-mu,mu,0,0),&tmp(0,0));
00275       y0=x-mu;
00276       y1=y0-nu;
00277       y2=y1+mu;
00278       _sse_mulABC_set_333(&U(y1,nu,0,0),&U(y0,mu,0,0),&b1(0,0));
00279       _sse_mulAHBC_set_333(&U(y1,mu,0,0),&b1(0,0),&b2(0,0));
00280       _sse_mulAHBC_set_333(&U(y2,mu,0,0),&b2(0,0),&b3(0,0));
00281       _sse_mulAHBC_add_333(&U(y2+mu,nu,0,0),&b3(0,0),&tmp(0,0));
00282       y0=x+mu;
00283       y1=y0+nu;
00284       y2=y1-mu;
00285       _sse_mulABC_set_333(&U(y0,mu,0,0),&U(y0+mu,nu,0,0),&b1(0,0));
00286       _sse_mulABHC_set_333(&b1(0,0),&U(y1,mu,0,0),&b2(0,0));
00287       _sse_mulABHC_set_333(&b2(0,0),&U(y2,mu,0,0),&b3(0,0));
00288       _sse_mulABHC_add_333(&b3(0,0),&U(x,nu,0,0),&tmp(0,0));
00289       y0=((x+mu)+mu)-nu;
00290       y1=y0-mu;
00291       y2=y1-mu;
00292       _sse_mulABHC_set_333(&U(x+mu,mu,0,0),&U(y0,nu,0,0),&b1(0,0));
00293       _sse_mulABHC_set_333(&b1(0,0),&U(y1,mu,0,0),&b2(0,0));
00294       _sse_mulABHC_set_333(&b2(0,0),&U(y2,mu,0,0),&b3(0,0));
00295       _sse_mulABC_add_333(&b3(0,0),&U(y2,nu,0,0),&tmp(0,0));
00296     }
00297     return tmp;
00298   }
00299 
00300   // //////////////////////////////////////////////////////
00301   // this is slow but should make the chair correcly ...
00302   // see: hep-lat/0712010
00303   // //////////////////////////////////////////////////////
00304 
00305   static mdp_matrix chair_H(gauge_field &U, site x, int mu) {
00306     int nc=3;
00307     int ndim=U.ndim;
00308     int nu, rho;
00309     mdp_matrix tmp(nc, nc);
00310     mdp_matrix b1(nc, nc);
00311     mdp_matrix b2(nc, nc);
00312     mdp_matrix b3(nc, nc);
00313     site y1(U.lattice());
00314     site y2(U.lattice());
00315     site y3(U.lattice());
00316     site y4(U.lattice());
00317     site y5(U.lattice());
00318     tmp=0;
00319     for(nu=0; nu<ndim; nu++) if(nu!=mu)
00320       for(rho=0; rho<ndim; rho++) if((rho!=nu) && (rho!=mu)) { 
00321         y1=x+mu;
00322         y2=y1+nu;
00323         y3=y2+rho;
00324         y4=y3-mu;
00325         y5=y4-nu;
00326         _sse_mulABC_set_333(&U(y1,nu,0,0),&U(y2,rho,0,0),&b1(0,0));
00327         _sse_mulABHC_set_333(&b1(0,0),&U(y4,mu,0,0),&b2(0,0));
00328         _sse_mulABHC_set_333(&b2(0,0),&U(y5,nu,0,0),&b3(0,0));
00329         _sse_mulABHC_add_333(&b3(0,0),&U(x,rho,0,0),&tmp(0,0));
00330         y1=x+mu;
00331         y2=y1-nu;
00332         y3=y2+rho;
00333         y4=y3-mu;
00334         y5=y4+nu;
00335         _sse_mulAHBC_set_333(&U(y2,nu,0,0),&U(y2,rho,0,0),&b1(0,0));
00336         _sse_mulABHC_set_333(&b1(0,0),&U(y4,mu,0,0),&b2(0,0));
00337         _sse_mulABC_set_333(&b2(0,0),&U(y4,nu,0,0),&b3(0,0));
00338         _sse_mulABHC_add_333(&b3(0,0),&U(x,rho,0,0),&tmp(0,0));
00339         y1=x+mu;
00340         y2=y1+nu;
00341         y3=y2-rho;
00342         y4=y3-mu;
00343         y5=y4-nu;
00344         _sse_mulABHC_set_333(&U(y1,nu,0,0),&U(y3,rho,0,0),&b1(0,0));
00345         _sse_mulABHC_set_333(&b1(0,0),&U(y4,mu,0,0),&b2(0,0));
00346         _sse_mulABHC_set_333(&b2(0,0),&U(y5,nu,0,0),&b3(0,0));
00347         _sse_mulABC_add_333(&b3(0,0),&U(y5,rho,0,0),&tmp(0,0));
00348         y1=x+mu;
00349         y2=y1-nu;
00350         y3=y2-rho;
00351         y4=y3-mu;
00352         y5=y4+nu;
00353         _sse_mulABC_set_333(&U(y4,mu,0,0),&U(y3,rho,0,0),&b1(0,0));
00354         _sse_mulABC_set_333(&b1(0,0),&U(y2,nu,0,0),&b2(0,0));
00355         _sse_mulAHBC_set_333(&b2(0,0),&U(y4,nu,0,0),&b3(0,0));
00356         _sse_mulABC_add_333(&b3(0,0),&U(y5,rho,0,0),&tmp(0,0));
00357       }
00358     return(tmp);
00359   }
00360   
00361   
00362   static mdp_matrix twisted_rectangle_H(gauge_field &U, site x, int mu) {
00363     int nu;
00364     int nc=3;
00365     site y1(U.lattice());
00366     site y2(U.lattice());
00367     site y3(U.lattice());
00368     site y4(U.lattice());
00369     site y5(U.lattice());
00370     mdp_matrix tmp(nc,nc),b1(nc,nc),b2(nc,nc);
00371     tmp=0;
00372     b1=0;
00373     b2=0;
00374     for(nu=0; nu<U.ndim; nu++)
00375       if(nu!=mu) {
00376         
00377         // type (a) staples /////////////////////////////////////////////
00378         
00379         y1=x+mu;
00380         y2=y1+mu;
00381         y4=x+nu;
00382         y3=y4+mu;
00383         _sse_mulABC_set_333(&U(y1,nu,0,0),&U(y3,mu,0,0),&b1(0,0));
00384         _sse_mulABHC_set_333(&b1(0,0), &U(y2,nu,0,0),&b2(0,0));
00385         _sse_mulABHC_set_333(&b2(0,0), &U(y1,mu,0,0),&b1(0,0));
00386         _sse_mulABC_set_333(&b1(0,0), &U(y1,nu,0,0),&b2(0,0));
00387         _sse_mulABHC_set_333(&b2(0,0), &U(y4,mu,0,0),&b1(0,0));
00388         _sse_mulABHC_add_333(&b1(0,0), &U(x,nu,0,0),&tmp(0,0));
00389         y2=y1-nu;
00390         y3=y2+mu;
00391         y4=y2-mu;
00392         _sse_mulAHBC_set_333(&U(y2,nu,0,0),&U(y2,mu,0,0),&b1(0,0));
00393         _sse_mulABC_set_333(&b1(0,0), &U(y3,nu,0,0),&b2(0,0));
00394         _sse_mulABHC_set_333(&b2(0,0), &U(y1,mu,0,0),&b1(0,0));
00395         _sse_mulABHC_set_333(&b1(0,0), &U(y2,nu,0,0),&b2(0,0));
00396         _sse_mulABHC_set_333(&b2(0,0), &U(y4,mu,0,0),&b1(0,0));
00397         _sse_mulABC_add_333(&b1(0,0), &U(y4,nu,0,0),&tmp(0,0));
00398         
00399         //   type (b) staples  //////////////////////////////////////////          
00400         
00401         y2=y1+nu;
00402         y3=x+nu;
00403         y4=y3+nu;
00404         _sse_mulABHC_set_333(&U(y1,nu,0,0),&U(y3,mu,0,0),&b1(0,0));
00405         _sse_mulABC_set_333(&b1(0,0), &U(y3,nu,0,0),&b2(0,0));
00406         _sse_mulABC_set_333(&b2(0,0), &U(y4,mu,0,0),&b1(0,0));
00407         _sse_mulABHC_set_333(&b1(0,0), &U(y2,nu,0,0),&b2(0,0));
00408         _sse_mulABHC_set_333(&b2(0,0), &U(y3,mu,0,0),&b1(0,0));
00409         _sse_mulABHC_add_333(&b1(0,0), &U(x,nu,0,0),&tmp(0,0));
00410         y2=x-nu;
00411         y3=y2-nu;
00412         y4=y3+mu;
00413         y5=y4+nu;
00414         _sse_mulABC_set_333(&U(y2,mu,0,0),&U(y5,nu,0,0),&b1(0,0));
00415         _sse_mulABC_set_333(&U(y3,nu,0,0),&b1(0,0),&b2(0,0));
00416         _sse_mulAHBC_set_333(&b2(0,0), &U(y3,mu,0,0),&b1(0,0));
00417         _sse_mulABC_set_333(&b1(0,0), &U(y4,nu,0,0),&b2(0,0));
00418         _sse_mulABHC_set_333(&b2(0,0), &U(y2,mu,0,0),&b1(0,0));
00419         _sse_mulABC_add_333(&b1(0,0), &U(y2,nu,0,0),&tmp(0,0)); 
00420         
00421         // type (c) staples /////////////////////////////////////////////
00422         
00423         y2=x+nu;
00424         y3=y2-mu;
00425         y4=y3-nu;
00426         _sse_mulABHC_set_333(&U(y1,nu,0,0),&U(y2,mu,0,0),&b1(0,0));
00427         _sse_mulABHC_set_333(&b1(0,0), &U(x,nu,0,0),&b2(0,0));
00428         _sse_mulABHC_set_333(&b2(0,0), &U(y4,mu,0,0),&b1(0,0));
00429         _sse_mulABC_set_333(&b1(0,0), &U(y4,nu,0,0),&b2(0,0));
00430         _sse_mulABC_set_333(&b2(0,0), &U(y3,mu,0,0),&b1(0,0));
00431         _sse_mulABHC_add_333(&b1(0,0), &U(x,nu,0,0),&tmp(0,0)); 
00432         y2=x-nu;
00433         y3=y2-mu;
00434         y4=y3+nu;
00435         y5=y2+mu;
00436         _sse_mulABC_set_333(&U(y2,mu,0,0),&U(y5,nu,0,0),&b1(0,0));
00437         _sse_mulAHBC_set_333(&b1(0,0), &U(y2,nu,0,0),&b2(0,0));
00438         _sse_mulABHC_set_333(&b2(0,0), &U(y4,mu,0,0),&b1(0,0));
00439         _sse_mulABHC_set_333(&b1(0,0), &U(y3,nu,0,0),&b2(0,0));
00440         _sse_mulABC_set_333(&b2(0,0), &U(y3,mu,0,0),&b1(0,0));
00441         _sse_mulABC_add_333(&b1(0,0), &U(y2,nu,0,0),&tmp(0,0));
00442         
00443         // type (d) staples /////////////////////////////////////////////
00444         
00445         y2=y1-nu;
00446         y3=y2-mu;
00447         y4=x+nu;
00448         _sse_mulABHC_set_333(&U(y1,nu,0,0),&U(y4,mu,0,0),&b1(0,0));
00449         _sse_mulABHC_set_333(&b1(0,0), &U(x,nu,0,0),&b2(0,0));
00450         _sse_mulABC_set_333(&b2(0,0), &U(x,mu,0,0),&b1(0,0));
00451         _sse_mulABHC_set_333(&b1(0,0), &U(y2,nu,0,0),&b2(0,0));
00452         _sse_mulABHC_set_333(&b2(0,0), &U(y3,mu,0,0),&b1(0,0));
00453         _sse_mulABC_add_333(&b1(0,0), &U(y3,nu,0,0),&tmp(0,0));
00454         
00455         // To include the next set double counts !!!    //////////////
00456         
00457       }
00458     return tmp;
00459   }
00460 
00461 #endif
00462 
00463 
00464 
00465   
00466 
00467   
00468   // ////////////////////////////////////////////////////////////////////
00469   // new_heatbath uses an improved gauge action!
00470   // both isotropic (param.zeta=1) and anisotropic
00471   // ////////////////////////////////////////////////////////////////////
00472 
00473   enum { iGauge_min=4 };
00474 
00475   static int strange_mapping(site &x) {
00476     int mu, type=0;
00477     for(mu=0; mu<x.lattice().ndim; mu++) type+=(int) pow((float) iGauge_min,mu)*(x(mu) % iGauge_min);
00478     return type;
00479   }
00480   
00481  public:
00482   static gauge_stats heatbath(gauge_field &U,
00483                               coefficients &coeff,
00484                               int n_iter=1,
00485                               string model="MILC") { 
00486     
00487     begin_function("ImprovedGaugeAction__heatbath");
00488 
00489     gauge_stats stats;
00490     mdp_real beta, zeta, u_s, u_t;
00491 
00492     if(coeff.has_key("beta")) beta=coeff["beta"]; else error("beta undeclared");
00493     if(coeff.has_key("zeta")) zeta=coeff["zeta"]; else zeta=1;
00494     if(coeff.has_key("u_t"))  u_t=coeff["u_t"];   else u_t=1;
00495     if(coeff.has_key("u_s"))  u_s=coeff["u_s"];   else u_s=1;
00496     
00497     // if(Nproc!=1)  error("improved_heatbath() does not work in parallel!");
00498 
00499     if(U.ndim!=4) error("fermiqcd_gauge_algorithms/improved_heatbath(): ndim!=4 (use heatbath instead)");
00500 
00501     /*
00502       if((U.lattice().size(0) % 3!=0) || (U.lattice().size(1) % 3!=0) ||
00503       (U.lattice().size(2) % 3!=0) || (U.lattice().size(3) % 3!=0))
00504       error("lattice is not divisible by 3");
00505     */
00506 
00507     if(U.nc==1)   error("fermiqcd_gauge_algorithms/improved_heatbath(): U(1)? (use metropolis instead)");
00508     int nc=U.nc;
00509     int ndim=U.ndim;
00510     int i,j,k,iter,mu,type;
00511     mdp_matrix M;
00512     site x(U.lattice());
00513     double time=mpi.time();
00514     mdp_complex a[4],tmpUik;
00515     mdp_real alpha_s;
00516     mdp_real c_tp=0, c_tr=0, c_sp=0, c_sr=0, c_p=0, c_r=0, c_c=0;
00517     
00518     if(model=="Morningstar") {
00519       
00520       c_tp=4.0*zeta/(3.0*pow((double)u_s*u_t,(double)2.0));
00521       c_tr=-1.0*zeta/(12.0*pow((double)u_s,(double)4.0)*pow((double)u_t,(double)2.0));
00522       c_sp=5.0/(3.0*zeta*pow((double)u_s,(double)4.0));
00523       c_sr=-1.0/(12.0*zeta*pow((double)u_s,(double)6.0));
00524       
00525       c_p=1.0*pow((double)u_s,(double)-4.0);
00526       c_r=-0.05*pow((double)u_s,(double)-6.0);
00527       c_c=0;
00528       
00529     } else if(model=="MILC") {
00530       
00531       if(zeta!=1) 
00532         error("fermiqcd_gauge_algorithms/improved_heatbath: zeta!=1");
00533       
00534       alpha_s=-4.0*log(u_s)/3.0684;
00535       c_p=1.0;
00536       c_r=-0.05*pow((double)u_s,(double)-2.0)*(1.0+0.4805*alpha_s);;
00537       c_c=-1.00*pow((double)u_s,(double)-2.0)*(0.03325*alpha_s);;
00538       
00539     } else {
00540       mdp << "Using default non-improved action" << '\n';
00541       stats=WilsonGaugeAction::heatbath(U,coeff,n_iter);
00542       end_function("ImprovedGaugeAction__heatbath");
00543       return stats;
00544     }
00545     
00546     mdp << coeff;
00547     
00548     for(iter=0; iter<n_iter; iter++)
00549       for(type=0; type<(int) pow((float) iGauge_min, U.ndim); type++) {
00550         forallsites(x) {
00551           if(strange_mapping(x)==type) {
00552             for(mu=0; mu<ndim; mu++) 
00553               for(i=0; i<nc-1; i++)
00554                 for(j=i+1; j<nc; j++) { 
00555                   if(zeta!=1) {
00556                     // anisotropic San Diego action
00557                     if(mu==0) M=U(x,mu)*(c_tp*staple_0i_H(U,x,0)+
00558                                          c_tr*rectangles_0i_H(U,x,0));
00559                     else      M=U(x,mu)*(c_sp*staple_ij_H(U,x,mu)+
00560                                          c_tp*staple_0i_H(U,x,mu)+
00561                                          c_sr*rectangles_ij_H(U,x,mu)+
00562                                          c_tr*rectangles_0i_H(U,x,mu));
00563                   } else if(c_c==0) {
00564                     // isotropic San Diego action
00565                     M=U(x,mu)*(c_p*staple_H(U,x,mu)+
00566                                c_r*rectangles_ij_H(U,x,mu,0));
00567                   } else {
00568                     // fully O(a^2) improved isotropic MILC action
00569                     M=U(x,mu)*(c_p*staple_H(U,x,mu)+
00570                                c_r*rectangles_ij_H(U,x,mu,0)+
00571                                c_c*chair_H(U,x,mu));
00572                   }
00573                   a[0]=M(i,i); 
00574                   a[1]=M(i,j);
00575                   a[2]=M(j,i);
00576                   a[3]=M(j,j);
00577                   heatbath_SU2(U.lattice().random(x),beta/U.nc,a);
00578                   for(k=0; k<U.nc; k++) {
00579                     tmpUik=a[0]*U(x,mu,i,k)+a[1]*U(x,mu,j,k);
00580                     U(x,mu,j,k)=a[2]*U(x,mu,i,k)+a[3]*U(x,mu,j,k);
00581                     U(x,mu,i,k)=tmpUik;
00582                   }
00583                 }
00584           }
00585         }
00586         U.update();
00587       }
00588     mdp << "\t<stats>\n\t\t<time>" << mpi.time()-time << "</time>\n\t</stats>\n";
00589     end_function("ImprovedGaugeAction__heatbath");
00590     return stats;
00591   }
00592 };
00593