fermiqcd_staggered_algorithms.h Source File

00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 mdp_matrix Omega4x4(mdp_site x) {
00014   mdp_matrix M(4,4);
00015   M=1;
00016   if(x(3) % 2==1) M=Gamma[3]*M;
00017   if(x(2) % 2==1) M=Gamma[2]*M;
00018   if(x(1) % 2==1) M=Gamma[1]*M;
00019   if(x(0) % 2==1) M=Gamma[0]*M;
00020   return M;
00021 }
00022 
00024 void (*default_staggered_action)(staggered_field &,
00025                                  staggered_field &,
00026                                  gauge_field &, 
00027                                  coefficients &, int) = StaggeredAsqtadActionFast::mul_Q;
00028 
00030 void mul_Q(staggered_field &psi_out,
00031            staggered_field &psi_in,
00032            gauge_field &U, 
00033            coefficients &coeff,
00034            int parity=EVENODD) {
00035   (*default_staggered_action)(psi_out, psi_in, U, coeff, parity);
00036 }
00037 
00038 // ////////////////////////////////////////////////
00039 // choice of the default inversion method
00040 // ////////////////////////////////////////////////
00041 
00043 inversion_stats (*default_staggered_inverter)(staggered_field &, 
00044                                               staggered_field &,
00045                                               gauge_field &, 
00046                                               coefficients &,
00047                                               mdp_real, mdp_real,int) 
00048   =BiCGStab::inverter<staggered_field,gauge_field>;
00049 
00051 inversion_stats mul_invQ(staggered_field &psi_out, 
00052                          staggered_field &psi_in, 
00053                          gauge_field &U, 
00054                          coefficients &coeff, 
00055                          mdp_real absolute_precision=staggered_inversion_precision,
00056                          mdp_real relative_precision=0, int max_steps=2000) {
00057   return (*default_staggered_inverter)(psi_out, psi_in, U, coeff, absolute_precision, relative_precision, max_steps);
00058 }
00059 
00060 // ////////////////////////////////////////////////////////////
00061 // ////////////////////////////////////////////////////////////
00062 // Lepage improved fat links for improved staggered action
00063 // ////////////////////////////////////////////////////////////
00064 // ////////////////////////////////////////////////////////////
00065 
00069 mdp_array<mdp_real,1> lepage_coefficients(mdp_real plaquette, char type[]) {
00070   begin_function("lepage_coefficients");
00071 
00072   mdp_real u0=pow((double)plaquette,(double)0.25);
00073   mdp_array<mdp_real,1> c(6); 
00074 
00075   c[0]=c[1]=c[2]=c[3]=c[4]=c[5]=0.0;
00076   if(strcmp(type,"None")==0) {
00077     c[0]=1;
00078   }
00079   if(strcmp(type,"Staple+Naik")==0) {
00080     c[0]=9.0/(8*4);
00081     c[1]=9.0/(8*8)*pow(u0,-2);
00082     c[2]=c[3]=c[4]=0;
00083     c[5]=-9.0/(8*27);
00084   }
00085   if(strcmp(type,"Fat3")==0) {
00086     error("Not implemented");
00087   }
00088   if(strcmp(type,"Fat5")==0) {
00089     c[0]=1.0/7;
00090     c[1]=1.0/(7*2)*pow(u0,-2);
00091     c[2]=1.0/(7*8)*pow(u0,-4);
00092     c[3]=c[4]=0;
00093     c[5]=0;
00094   }
00095   if(strcmp(type,"Fat7")==0) {
00096     c[0]=1.0/8;
00097     c[1]=1.0/(8*2)*pow(u0,-2);
00098     c[2]=1.0/(8*8)*pow(u0,-4);
00099     c[3]=1.0/(8*48)*pow(u0,-6);
00100     c[4]=0;
00101     c[5]=0;
00102   }
00103   if(strcmp(type,"Full")==0) {
00104     c[0]=5.0/8.0;
00105     c[1]=1.0/16.0*pow(u0,-2); 
00106     c[2]=1.0/64.0*pow(u0,-4);
00107     c[3]=1.0/(8.0*48.0)*pow(u0,-6);
00108     c[4]=-1.0/16.0*pow(u0,-4);
00109     c[5]=-1.0/24.0*pow(u0,-2);
00110   }
00111   
00112   mdp << "\tImprovement type=" << type << '\n';
00113   mdp << "Coefficients:\n";
00114   mdp << "u_0=" << u0 << '\n';
00115   mdp << "link=" << c[0] << '\n';
00116   mdp << "3staple=" << c[1] << '\n';
00117   mdp << "5staple=" << c[2] << '\n';
00118   mdp << "7staple=" << c[3] << '\n';
00119   mdp << "lepage=" << c[4] << '\n';
00120   mdp << "naik=" << c[5] << '\n';
00121 
00122   end_function("lepage_coefficients");
00123 
00124   return c;
00125 }
00126 
00127 
00128 
00129 // ///////////////////////////////////////////////
00130 // Note we cannot use Nmdp_matrixField here
00131 // because it does not implement the twist
00132 // we use vectors of gauge_fields instead!
00133 // ///////////////////////////////////////////////
00134 
00155 void lepage_improved_links(gauge_field &V, 
00156                            gauge_field &U, 
00157                            mdp_array<mdp_real,1> c, 
00158                            int project=FALSE) {
00159 
00160   begin_function("lepage_improved_links");
00161 
00162   if(U.ndim!=4) error("fermiqcd_staggered_auxiliary_functions/lepage_improved_links: ndim!=4 (contact <mdp@fnal.gov>)");
00163 
00164   // this is a list of permitations of 0,1,2,3
00165 
00166   
00167   const int epsilon[24][5] = {{0, 1, 2, 3, 0}, {0, 1, 3, 2, 3},
00168                               {0, 2, 1, 3, 1}, {0, 2, 3, 1, 4},
00169                               {0, 3, 1, 2, 2}, {0, 3, 2, 1, 5},
00170                               {1, 0, 2, 3, 0}, {1, 0, 3, 2, 3},
00171                               {1, 2, 0, 3, 1}, {1, 2, 3, 0, 4},
00172                               {1, 3, 0, 2, 2}, {1, 3, 2, 0, 5},
00173                               {2, 0, 1, 3, 0}, {2, 0, 3, 1, 3},
00174                               {2, 1, 0, 3, 1}, {2, 1, 3, 0, 4},
00175                               {2, 3, 0, 1, 2}, {2, 3, 1, 0, 5},
00176                               {3, 0, 1, 2, 0}, {3, 0, 2, 1, 3},
00177                               {3, 1, 0, 2, 1}, {3, 1, 2, 0, 4},
00178                               {3, 2, 0, 1, 2}, {3, 2, 1, 0, 5}}; 
00179   // mu, nu, rho, sig, im2 
00180 
00181   
00182   int nc=U.nc;
00183   int ndim=U.ndim;
00184   int mu,nu,rho, idx, imn, im2, i,j;
00185   site x(U.lattice());
00186   site y(U.lattice());
00187   mdp_matrix b1(nc,nc), b2(nc,nc);
00188 
00189   mdp << "Allocating temporary vectors...";
00190 
00191   gauge_field *Delta1=new gauge_field[(ndim-1)];
00192   gauge_field *Delta2=new gauge_field[(ndim-1)*(ndim-2)];
00193 
00194   mdp << "done.\n";
00195 
00196   for(imn=0; imn<(ndim-1); imn++)
00197     Delta1[imn].allocate_gauge_field(U.lattice(),nc);
00198   for(im2=0; im2<(ndim-1)*(ndim-2); im2++)
00199     Delta2[im2].allocate_gauge_field(U.lattice(),nc);
00200 
00201   // ////////////////////////////////////////////////////////////////////
00202   // These coefficints are in agreement with Irginos, Toussaint and Sugar
00203   // hep-lat/9903032
00204   // equivalent but easier than the Lepage version
00205   // ////////////////////////////////////////////////////////////////////
00206 
00207   
00208   mdp << "Computing Fat Links for Lepage improved action:\n";
00209 
00210   // LINK
00211   mdp << "link...\n";
00212 
00213   forallsites(x)
00214     for(mu=0; mu<ndim; mu++)
00215       V(x,mu)=c[0]*U(x,mu);
00216 
00217   // 3 STAPLE
00218   mdp << "3staple...\n";
00219 
00220   forallsites(x)
00221     for(mu=0; mu<ndim; mu++)
00222       for(nu=0; nu<ndim; nu++) if(nu!=mu) {
00223         imn=(nu<mu)?nu:(nu-1);
00224         Delta1[imn](x,mu) =U(x,nu)*U(x+nu,mu)*hermitian(U(x+mu,nu));
00225         y=x-nu;
00226         Delta1[imn](x,mu)+=hermitian(U(y,nu))*U(y,mu)*U(y+mu,nu);
00227         V(x,mu)+=c[1]*Delta1[imn](x,mu);
00228       }
00229 
00230   
00231   for(imn=0; imn<(ndim-1); imn++)
00232     Delta1[imn].update();  
00233 
00234   // 5 STAPLE
00235   mdp << "5staple...\n";
00236 
00237   forallsites(x) 
00238     for(idx=0; idx<24; idx++) {
00239       mu =epsilon[idx][0];
00240       nu =epsilon[idx][1];
00241       rho=epsilon[idx][2];
00242       im2=epsilon[idx][4];
00243       imn=(nu<mu)?nu:(nu-1);
00244       Delta2[im2](x,mu) =U(x,rho)*Delta1[imn](x+rho,mu)*hermitian(U(x+mu,rho));
00245       y=x-rho;
00246       Delta2[im2](x,mu)+=hermitian(U(y,rho))*Delta1[imn](y,mu)*U(y+mu,rho);
00247       V(x,mu)+=c[2]*Delta2[im2](x,mu);
00248     }
00249   for(im2=0; im2<(ndim-1)*(ndim-2); im2++)
00250     Delta2[im2].update();
00251   
00252   // 7 STAPLE
00253   mdp << "7staple...\n";
00254 
00255   if(c[3]!=0) 
00256     forallsites(x) 
00257       for(idx=0; idx<24; idx++) {
00258         mu =epsilon[idx][0];
00259         rho=epsilon[idx][3];
00260         im2=epsilon[idx][4];
00261         b2=U(x,rho)*Delta2[im2](x+rho,mu)*hermitian(U(x+mu,rho));
00262         y=x-rho;
00263         b2+=hermitian(U(y,rho))*Delta2[im2](y,mu)*U(y+mu,rho);
00264         V(x,mu)+=c[3]*b2;
00265       }
00266   
00267   // LEPAGE TERM UP
00268   mdp << "lepage...\n";
00269 
00270   if(c[4]!=0) {
00271     forallsites(x)
00272       for(mu=0; mu<ndim; mu++)
00273         for(nu=0; nu<ndim; nu++) if(nu!=mu) {
00274           imn=(nu<mu)?nu:(nu-1);
00275           Delta1[imn](x,mu)=U(x,nu)*U(x+nu,mu)*hermitian(U(x+mu,nu));
00276         }
00277     for(imn=0; imn<(ndim-1); imn++)
00278       Delta1[imn].update();
00279 
00280     forallsites(x)
00281       for(mu=0; mu<ndim; mu++)
00282         for(nu=0; nu<ndim; nu++) if(nu!=mu) {
00283           imn=(nu<mu)?nu:(nu-1);
00284           b2=U(x,nu)*Delta1[imn](x+nu,mu)*hermitian(U(x+mu,nu));
00285           V(x,mu)+=c[4]*b2;
00286         }
00287     
00288     // LEPAGE TERM DOWN
00289     forallsites(x)
00290       for(mu=0; mu<ndim; mu++)
00291         for(nu=0; nu<ndim; nu++) if(nu!=mu) {
00292           imn=(nu<mu)?nu:(nu-1);
00293           y=x-nu;
00294           Delta1[imn](x,mu)=hermitian(U(y,nu))*U(y,mu)*U(y+mu,nu);
00295         }
00296     for(imn=0; imn<(ndim-1); imn++)
00297       Delta1[imn].update();
00298     
00299     forallsites(x)
00300       for(mu=0; mu<ndim; mu++)
00301         for(nu=0; nu<ndim; nu++) if(nu!=mu) {
00302           imn=(nu<mu)?nu:(nu-1);
00303           y=x-nu;
00304           b2=hermitian(U(y,nu))*Delta1[imn](y,mu)*U(y+mu,nu);
00305           V(x,mu)+=c[4]*b2;
00306         }
00307   }
00308   if(project==TRUE) {
00309     mdp << "projecting to SU(n)...\n";
00310 
00311     forallsites(x)
00312       for(mu=0; mu<4; mu++)
00313         V(x,mu)=project_SU(V(x,mu));
00314   }
00315   V.update();
00316 
00317   if(c[5]!=0) {
00318     mdp << "naik...\n";
00319     compute_long_links(V,U,3);
00320     mdp << "normalizing naik...\n";
00321 
00322     forallsitesandcopies(x)
00323       for(mu=0; mu<U.ndim; mu++)
00324         for(i=0; i<U.nc; i++)
00325           for(j=0; j<U.nc; j++)
00326             V.long_links(x,mu,i,j)*=c[5];
00327   }
00328 
00329   cout << "Freeing temporary vectors...";
00330 
00331   for(imn=0; imn<(ndim-1); imn++)
00332     Delta1[imn].deallocate_memory();
00333   for(im2=0; im2<(ndim-1)*(ndim-2); im2++)
00334     Delta2[im2].deallocate_memory();
00335 
00336   cout << "done.\n";
00337 
00338   end_function("lepage_improved_links");
00339 }
00340 
00341 void staggered_rephase(gauge_field &U, staggered_field &chi) {
00342   
00343   begin_function("staggered_rephase");
00344   
00345   site x(U.lattice());
00346   int  mu,i,j;
00347   forallsites(x) {
00348     for(mu=0; mu<U.ndim; mu++) {
00349       for(i=0; i<U.nc; i++)
00350         for(j=0; j<U.nc; j++) 
00351           U(x,mu,i,j)=U(x,mu,i,j)*chi.eta(x,mu);
00352     }
00353     // this takes car of antiperiodic bc
00354     if(x(0)==U.lattice().size(0)-1) {
00355       for(i=0; i<U.nc; i++)
00356         for(j=0; j<U.nc; j++) 
00357           U(x,0,i,j)*=-1;
00358     }
00359   }
00360 
00361   end_function("staggered_rephase");
00362 
00363 }
00364