fermiqcd_gauge_routines.h Source File

00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 // ///////////////////////////
00014 // If use class mdp_matrix
00015 // ///////////////////////////
00016 inline mdp_matrix staple(gauge_field &U, register site x, 
00017                      int mu, int s1, int nu) {
00018   mdp_matrix tmp(U.nc,U.nc);
00019     if(s1==+1) {
00020       tmp=U(x,nu)*U(x+nu,mu)*hermitian(U(x+mu,nu));
00021     } else {
00022       site y(U.lattice());
00023       y=x-nu;
00024       tmp=hermitian(U(y,nu))*U(y,mu)*U(y+mu,nu);
00025     }
00026     return tmp;
00027 }
00028 inline mdp_matrix staple(gauge_field &U, register site x, int mu) {
00029   mdp_matrix tmp(U.nc,U.nc);
00030   site y(U.lattice());
00031   int nu;
00032   tmp=0;
00033     for(nu=0; nu<U.ndim; nu++) if(nu!=mu) {
00034       tmp+=U(x,nu)*U(x+nu,mu)*hermitian(U(x+mu,nu));
00035       y=x-nu;
00036       tmp+=hermitian(U(y,nu))*U(y,mu)*U(y+mu,nu);
00037     }
00038     return tmp;
00039   }
00040 
00041 inline mdp_matrix staple_H(gauge_field &U, register site x, 
00042                        int mu, int s1, int nu) {
00043   mdp_matrix tmp(U.nc,U.nc);
00044   if(s1==+1) {
00045     tmp=U(x+mu,nu)*hermitian(U(x+nu,mu))*hermitian(U(x,nu));
00046   } else {
00047       site y(U.lattice());
00048       y=x-nu;
00049       tmp=hermitian(U(y+mu,nu))*hermitian(U(y,mu))*U(y,nu);
00050   }
00051   return tmp;
00052 }
00053 inline mdp_matrix staple_H(gauge_field &U, register site x, int mu) {
00054   mdp_matrix tmp(U.nc,U.nc);
00055   site y(U.lattice());
00056   int nu;
00057   tmp=0;
00058   for(nu=0; nu<U.ndim; nu++) if(nu!=mu) {
00059     tmp+=U(x+mu,nu)*hermitian(U(x+nu,mu))*hermitian(U(x,nu));
00060     y=x-nu;
00061     tmp+=hermitian(U(y+mu,nu))*hermitian(U(y,mu))*U(y,nu);
00062   }
00063   return tmp;
00064 }
00065 
00066 inline mdp_matrix staple_0i_H(gauge_field &U, register site x, int mu) {
00067   mdp_matrix tmp(U.nc,U.nc);
00068   site y(U.lattice());
00069   int nu;
00070   if(mu==0) {
00071     tmp=0;
00072     for(nu=1; nu<U.ndim; nu++) {
00073       tmp+=U(x+mu,nu)*hermitian(U(x+nu,mu))*hermitian(U(x,nu));
00074       y=x-nu;
00075       tmp+=hermitian(U(y+mu,nu))*hermitian(U(y,mu))*U(y,nu);
00076     }
00077   } else {
00078     nu=0;
00079     tmp=U(x+mu,nu)*hermitian(U(x+nu,mu))*hermitian(U(x,nu));
00080     y=x-nu;
00081     tmp+=hermitian(U(y+mu,nu))*hermitian(U(y,mu))*U(y,nu);
00082     }
00083   
00084   return tmp;
00085 }  
00086 
00087 inline mdp_matrix staple_ij_H(gauge_field &U, register site x, int mu) {
00088   mdp_matrix tmp(U.nc,U.nc);
00089   site y(U.lattice());
00090   int nu;
00091   tmp=0;
00092   if(mu!=0)
00093     for(nu=1; nu<U.ndim; nu++) if(nu!=mu) {
00094       tmp+=U(x+mu,nu)*hermitian(U(x+nu,mu))*hermitian(U(x,nu));
00095         y=x-nu;
00096         tmp+=hermitian(U(y+mu,nu))*hermitian(U(y,mu))*U(y,nu);
00097     }
00098   return tmp;
00099 }
00100 
00101 // ///////////////////////////
00102 // plaquette
00103 // ////////////////////////// 
00104 inline mdp_matrix plaquette(gauge_field &U, site x, int mu, int nu) {
00105   mdp_matrix tmp(U.nc,U.nc);
00106   tmp=U(x,mu)*U(x+mu,nu)*hermitian(U(x+nu,mu))*hermitian(U(x,nu));
00107   return tmp;
00108 }
00109 
00110 /* obsolete stuff 
00111 
00112 void fast_mul_AB_to_C(mdp_complex *a, mdp_complex *b, mdp_complex *c, 
00113                       int &ni) {
00114   register int i,j,k, ink, inj;
00115   for(i=0; i<ni; i++) {
00116     ink=i*ni; inj=i*ni;
00117     for(k=0; k<ni; k++) {
00118       c[ink+k]=a[inj]*b[k];
00119       for(j=1; j<ni; j++) c[ink+k]+=a[inj+j]*b[j*ni+k];
00120     }
00121   }
00122 }
00123 void fast_mul_ABH_to_C(mdp_complex *a, mdp_complex *b, mdp_complex *c, 
00124                        int &ni) {
00125   register int i,j,k, ink, inj;
00126   for(i=0; i<ni; i++) {
00127     ink=i*ni; inj=i*ni;
00128     for(k=0; k<ni; k++) {
00129       c[ink+k]=a[inj]*conj(b[k*ni]);
00130       for(j=1; j<ni; j++) c[ink+k]+=a[inj+j]*conj(b[k*ni+j]);
00131     }
00132   }
00133 }
00134 void fast_mul_AHB_to_C(mdp_complex *a, mdp_complex *b, mdp_complex *c, 
00135                        int &ni) {
00136   register int i,j,k, ink;
00137   for(i=0; i<ni; i++) {
00138     ink=i*ni;
00139     for(k=0; k<ni; k++) {
00140       c[ink+k]=conj(a[i])*b[k];
00141       for(j=1; j<ni; j++) c[ink+k]+=conj(a[j*ni+i])*b[j*ni+k];
00142     }
00143   }
00144 }
00145 void fast_mul_AHBH_to_C(mdp_complex *a, mdp_complex *b, mdp_complex *c, 
00146                        int &ni) {
00147   register int i,j,k, ink, inj;
00148   for(i=0; i<ni; i++) {
00149     ink=i*ni; inj=i*ni;
00150     for(k=0; k<ni; k++) {
00151       c[ink+k]=conj(a[i]*b[k*ni]);
00152       for(j=1; j<ni; j++) c[ink+k]+=conj(a[j*ni+i]*b[k*ni+j]);
00153     }
00154   }
00155 }
00156 
00157 
00158 void fast_mul_AB_addto_C(mdp_complex *a, mdp_complex *b, mdp_complex *c, 
00159                       int &ni) {
00160   register int i,j,k, ink, inj;
00161   for(i=0; i<ni; i++) {
00162     ink=i*ni; inj=i*ni;
00163     for(k=0; k<ni; k++) {
00164       c[ink+k]+=a[inj]*b[k];
00165       for(j=1; j<ni; j++) c[ink+k]+=a[inj+j]*b[j*ni+k];
00166     }
00167   }
00168 }
00169 void fast_mul_ABH_addto_C(mdp_complex *a, mdp_complex *b, mdp_complex *c, 
00170                        int &ni) {
00171   register int i,j,k, ink, inj;
00172   for(i=0; i<ni; i++) {
00173     ink=i*ni; inj=i*ni;
00174     for(k=0; k<ni; k++) {
00175       c[ink+k]+=a[inj]*conj(b[k*ni]);
00176       for(j=1; j<ni; j++) c[ink+k]+=a[inj+j]*conj(b[k*ni+j]);
00177     }
00178   }
00179 }
00180 void fast_mul_AHB_addto_C(mdp_complex *a, mdp_complex *b, mdp_complex *c, 
00181                        int &ni) {
00182   register int i,j,k, ink;
00183   for(i=0; i<ni; i++) {
00184     ink=i*ni;
00185     for(k=0; k<ni; k++) {
00186       c[ink+k]+=conj(a[i])*b[k];
00187       for(j=1; j<ni; j++) c[ink+k]+=conj(a[j*ni+i])*b[j*ni+k];
00188     }
00189   }
00190 }
00191 void fast_add_AB_to_C(mdp_complex *a, mdp_complex *b, mdp_complex *c, 
00192                       int &ni) {
00193   for(register int i=0; i<ni*ni; i++) c[i]=a[i]+b[i];
00194 }
00195 void fast_add_AB_addto_C(mdp_complex *a, mdp_complex *b, mdp_complex *c, 
00196                          int &ni) {
00197   for(register int i=0; i<ni*ni; i++) c[i]+=a[i]+b[i];
00198 }
00199 void fast_add_aAbB_to_C(mdp_complex c1, mdp_complex *a1, mdp_complex c2, mdp_complex *a2,
00200                         mdp_complex *c, int &ni) {
00201   for(register int i=0; i<ni*ni; i++) c[i]=c1*a1[i]+c2*a2[i];
00202 }
00203 void fast_add_aAbB_addto_C(mdp_complex c1, mdp_complex *a1, mdp_complex c2, mdp_complex *a2,
00204                         mdp_complex *c, int &ni) {
00205   for(register int i=0; i<ni*ni; i++) c[i]+=c1*a1[i]+c2*a2[i];
00206 }
00207 void fast_add_aAbB_to_C(mdp_complex c1, mdp_complex *a1, mdp_complex c2, mdp_complex *a2,
00208                         mdp_complex c3, mdp_complex *a3, mdp_complex *c, int &ni) {
00209   for(register int i=0; i<ni*ni; i++) c[i]=c1*a1[i]+c2*a2[i]+c3*a3[i];
00210 }
00211 void fast_add_aAbB_addto_C(mdp_complex c1, mdp_complex *a1, mdp_complex c2, mdp_complex *a2,
00212                         mdp_complex c3, mdp_complex *a3, mdp_complex *c, int &ni) {
00213   for(register int i=0; i<ni*ni; i++) c[i]+=c1*a1[i]+c2*a2[i]+c3*a3[i];
00214 }
00215 void fast_scalar_mul_AB_to_C(mdp_complex a, mdp_complex *b, mdp_complex *c, 
00216                              int &ni) {
00217   for(register int i=0; i<ni*ni; i++) c[i]=a*b[i];
00218 }
00219 void fast_scalar_mul_AB_addto_C(mdp_complex a, mdp_complex *b, mdp_complex *c, 
00220                                 int &ni) {
00221   for(register int i=0; i<ni*ni; i++) c[i]+=a*b[i];
00222 }
00223 
00224 
00225 
00226 
00227 void fast_mul_AB_subto_C(mdp_complex *a, mdp_complex *b, mdp_complex *c, 
00228                       int &ni) {
00229   register int i,j,k, ink, inj;
00230   for(i=0; i<ni; i++) {
00231     ink=i*ni; inj=i*ni;
00232     for(k=0; k<ni; k++) {
00233       c[ink+k]+=a[inj]*b[k];
00234       for(j=1; j<ni; j++) c[ink+k]-=a[inj+j]*b[j*ni+k];
00235     }
00236   }
00237 }
00238 void fast_mul_ABH_subto_C(mdp_complex *a, mdp_complex *b, mdp_complex *c, 
00239                        int &ni) {
00240   register int i,j,k, ink, inj;
00241   for(i=0; i<ni; i++) {
00242     ink=i*ni; inj=i*ni;
00243     for(k=0; k<ni; k++) {
00244       c[ink+k]+=a[inj]*conj(b[k*ni]);
00245       for(j=1; j<ni; j++) c[ink+k]-=a[inj+j]*conj(b[k*ni+j]);
00246     }
00247   }
00248 }
00249 void fast_mul_AHB_subto_C(mdp_complex *a, mdp_complex *b, mdp_complex *c, 
00250                        int &ni) {
00251   register int i,j,k, ink;
00252   for(i=0; i<ni; i++) {
00253     ink=i*ni;
00254     for(k=0; k<ni; k++) {
00255       c[ink+k]+=conj(a[i])*b[k];
00256       for(j=1; j<ni; j++) c[ink+k]-=conj(a[j*ni+i])*b[j*ni+k];
00257     }
00258   }
00259 }
00260 void fast_add_AB_subto_C(mdp_complex *a, mdp_complex *b, mdp_complex *c, 
00261                          int &ni) {
00262   for(register int i=0; i<ni*ni; i++) c[i]-=a[i]+b[i];
00263 }
00264 void fast_add_aAbB_subto_C(mdp_complex c1, mdp_complex *a1, mdp_complex c2, mdp_complex *a2,
00265                         mdp_complex *c, int &ni) {
00266   for(register int i=0; i<ni*ni; i++) c[i]-=c1*a1[i]+c2*a2[i];
00267 }
00268 void fast_add_aAbB_subto_C(mdp_complex c1, mdp_complex *a1, mdp_complex c2, mdp_complex *a2,
00269                         mdp_complex c3, mdp_complex *a3, mdp_complex *c, int &ni) {
00270   for(register int i=0; i<ni*ni; i++) c[i]-=c1*a1[i]+c2*a2[i]+c3*a3[i];
00271 }
00272 void fast_scalar_mul_AB_subto_C(mdp_complex a, mdp_complex *b, mdp_complex *c, 
00273                                 int &ni) {
00274   for(register int i=0; i<ni*ni; i++) c[i]-=a*b[i];
00275 }
00276 
00277 // ///////////////////////////
00278 // else use fast_mul functions
00279 // gains a factor 1.45 in time
00280   // ///////////////////////////
00281 
00282 inline mdp_matrix staple(gauge_field &U, register site x, 
00283                      int mu, int s1, int nu) {
00284   int nc=U.nc;
00285   mdp_matrix b1(nc,nc);
00286   site y(U.lattice());
00287   mdp_matrix tmp(nc,nc);
00288   if(s1==+1) {
00289     fast_mul_AB_to_C(&U(x,nu,0,0),&U(x+nu,mu,0,0),&b1(0,0),nc);
00290       fast_mul_ABH_to_C(&b1(0,0),&U(x+mu,nu,0,0),&tmp(0,0),nc);
00291   } else {
00292     y=x-nu;
00293     fast_mul_AB_to_C(&U(y,mu,0,0),&U(y+mu,nu,0,0),&b1(0,0),nc);
00294     fast_mul_AHB_to_C(&U(y,nu,0,0),&b1(0,0),&tmp(0,0),nc);
00295   }
00296     return tmp;
00297 }
00298 inline mdp_matrix staple(gauge_field &U, register site x, int mu) {
00299   int nc=U.nc;
00300     mdp_matrix b1(nc,nc);
00301     mdp_matrix tmp(nc,nc);
00302     site y(U.lattice());
00303     int nu;
00304     tmp=0;
00305     for(nu=0; nu<U.ndim; nu++) if(nu!=mu) {
00306       fast_mul_AB_to_C(&U(x,nu,0,0),&U(x+nu,mu,0,0),&b1(0,0),nc);
00307       fast_mul_ABH_addto_C(&b1(0,0),&U(x+mu,nu,0,0),&tmp(0,0),nc);
00308       y=x-nu;
00309       fast_mul_AB_to_C(&U(y,mu,0,0),&U(y+mu,nu,0,0),&b1(0,0),nc);
00310       fast_mul_AHB_addto_C(&U(y,nu,0,0),&b1(0,0),&tmp(0,0),nc);
00311     }
00312     return tmp;
00313 }
00314 
00315 inline mdp_matrix staple_H(gauge_field &U, register site x, 
00316                               int mu, int s1, int nu) {
00317   int nc=U.nc;
00318   mdp_matrix b1(nc,nc);
00319   mdp_matrix tmp(nc,nc);
00320   site y(U.lattice());
00321   if(s1==+1) {
00322     fast_mul_ABH_to_C(&U(x+mu,nu,0,0),&U(x+nu,mu,0,0),&b1(0,0),nc);
00323     fast_mul_ABH_to_C(&b1(0,0),&U(x,nu,0,0),&tmp(0,0),nc);
00324   } else {
00325     y=x-nu;
00326     fast_mul_AHB_to_C(&U(y,mu,0,0),&U(y,nu,0,0),&b1(0,0),nc);
00327     fast_mul_AHB_to_C(&U(y+mu,nu,0,0),&b1(0,0),&tmp(0,0),nc);
00328   }
00329   return tmp;
00330 }
00331 inline mdp_matrix staple_H(gauge_field &U, register site x, int mu) {
00332   int nc=U.nc;
00333   mdp_matrix b1(nc,nc);
00334   mdp_matrix tmp(nc,nc);
00335   site y(U.lattice());
00336   int nu;
00337   tmp=0;
00338   for(nu=0; nu<U.ndim; nu++) if(nu!=mu) {
00339     fast_mul_ABH_to_C(&U(x+mu,nu,0,0),&U(x+nu,mu,0,0),&b1(0,0),nc);
00340     fast_mul_ABH_addto_C(&b1(0,0),&U(x,nu,0,0),&tmp(0,0),nc);
00341     y=x-nu;
00342     fast_mul_AHB_to_C(&U(y,mu,0,0),&U(y,nu,0,0),&b1(0,0),nc);
00343     fast_mul_AHB_addto_C(&U(y+mu,nu,0,0),&b1(0,0),&tmp(0,0),nc);
00344   }
00345   return tmp;
00346 }
00347 
00348 inline mdp_matrix staple_0i_H(gauge_field &U, register site x, int mu) {
00349   int nc=U.nc;
00350   mdp_matrix b1(nc,nc);
00351   mdp_matrix tmp(U.nc,U.nc);
00352   site y(U.lattice());
00353   int nu;
00354   if(mu==0) {
00355     tmp=0;
00356     for(nu=1; nu<U.ndim; nu++) {
00357       fast_mul_ABH_to_C(&U(x+mu,nu,0,0),&U(x+nu,mu,0,0),&b1(0,0),nc);
00358       fast_mul_ABH_addto_C(&b1(0,0),&U(x,nu,0,0),&tmp(0,0),nc);
00359       y=x-nu;
00360       fast_mul_AHB_to_C(&U(y,mu,0,0),&U(y,nu,0,0),&b1(0,0),nc);
00361       fast_mul_AHB_addto_C(&U(y+mu,nu,0,0),&b1(0,0),&tmp(0,0),nc);
00362     }
00363   } else {
00364     nu=0;
00365     fast_mul_ABH_to_C(&U(x+mu,nu,0,0),&U(x+nu,mu,0,0),&b1(0,0),nc);
00366     fast_mul_ABH_to_C(&b1(0,0),&U(x,nu,0,0),&tmp(0,0),nc);
00367     y=x-nu;
00368     fast_mul_AHB_to_C(&U(y,mu,0,0),&U(y,nu,0,0),&b1(0,0),nc);
00369     fast_mul_AHB_addto_C(&U(y+mu,nu,0,0),&b1(0,0),&tmp(0,0),nc);
00370   }
00371   return tmp;
00372 }  
00373 
00374 inline mdp_matrix staple_ij_H(gauge_field &U, register site x, int mu) {
00375   int nc=U.nc;
00376   mdp_matrix b1(nc,nc);
00377   mdp_matrix tmp(U.nc,U.nc);
00378   site y(U.lattice());
00379   int nu;
00380   tmp=0;
00381   if(mu!=0)
00382     for(nu=1; nu<U.ndim; nu++) if(nu!=mu) {
00383       fast_mul_ABH_to_C(&U(x+mu,nu,0,0),&U(x+nu,mu,0,0),&b1(0,0),nc);
00384       fast_mul_ABH_addto_C(&b1(0,0),&U(x,nu,0,0),&tmp(0,0),nc);
00385       y=x-nu;
00386       fast_mul_AHB_to_C(&U(y,mu,0,0),&U(y,nu,0,0),&b1(0,0),nc);
00387       fast_mul_AHB_addto_C(&U(y+mu,nu,0,0),&b1(0,0),&tmp(0,0),nc);
00388     }
00389   return tmp;
00390 }
00391 
00392 // ///////////////////////////
00393 // plaquette
00394 // ////////////////////////// 
00395 inline mdp_matrix plaquette(gauge_field &U, site x, int mu, int nu) {
00396   int nc=U.nc;
00397   mdp_matrix b1(nc,nc);
00398   mdp_matrix b2(nc,nc);
00399   mdp_matrix tmp(U.nc,U.nc);
00400   fast_mul_AB_to_C(&U(x,mu,0,0),&U(x+mu,nu,0,0),&b1(0,0),nc);
00401   fast_mul_ABH_to_C(&b1(0,0),&U(x+nu,mu,0,0),&b2(0,0),nc);
00402   fast_mul_ABH_to_C(&b2(0,0),&U(x,nu,0,0),&tmp(0,0),nc);
00403   return tmp;
00404 }
00405 */
00406