fermiqcd: /Users/mdipierro/fermiqcd/development/Libraries/fermiqcd_staggered

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