21 #if defined USE_GROUP_SU3 
   22 #include "field_G_imp_SU3.inc" 
   23 #elif defined USE_GROUP_SU2 
   24 #include "field_G_imp_SU2.inc" 
   25 #elif defined USE_GROUP_SU_N 
   26 #include "field_G_imp_SU_N.inc" 
   50   int is = 
nvol() * i_thread / Nthread;
 
   51   int ns = 
nvol() * (i_thread + 1) / Nthread - is;
 
   53   for (
int mu = 0; mu < 
nex(); ++mu) {
 
   54     for (
int site = is; site < is + ns; ++site) {
 
   68   for (
int mu = 0; mu < 
m_Nex; ++mu) {
 
   69     for (
int site = 0; site < 
m_Nvol; ++site) {
 
   70       this->
mat(ut, site, mu);
 
   95   int is = 
nvol() * i_thread / Nthread;
 
   96   int ns = 
nvol() * (i_thread + 1) / Nthread - is;
 
   98   for (
int mu = 0; mu < 
nex(); ++mu) {
 
   99     for (
int site = is; site < is + ns; ++site) {
 
  113   assert(ex < W.
nex());
 
  114   assert(ex1 < U1.
nex());
 
  115   assert(ex2 < U2.
nex());
 
  124   int is = W.
nvol() * i_thread / Nthread;
 
  125   int ns = W.
nvol() * (i_thread + 1) / Nthread - is;
 
  127   double       *g  = W.
ptr(0, is, ex);
 
  128   const double *g1 = U1.
ptr(0, is, ex1);
 
  129   const double *g2 = U2.
ptr(0, is, ex2);
 
  131   const int Nc  = W.
nc();
 
  132   const int Nc2 = 2 * Nc;
 
  133   const int Ndf = Nc2 * Nc;
 
  135   for (
int site = 0; site < ns; ++site) {
 
  138     for (
int ic2 = 0; ic2 < Nc; ++ic2) {
 
  139       for (
int ic1 = 0; ic1 < Nc; ++ic1) {
 
  140         int jg2 = ic2 * 2 + ig;
 
  141         int jg1 = ic1 * Nc2 + ig;
 
  142         g[2 * ic2 + Nc2 * ic1 + ig]     = mult_Gnn_r(&g1[jg1], &g2[jg2], Nc);
 
  143         g[2 * ic2 + 1 + Nc2 * ic1 + ig] = mult_Gnn_i(&g1[jg1], &g2[jg2], Nc);
 
  155   assert(ex < W.
nex());
 
  156   assert(ex1 < U1.
nex());
 
  157   assert(ex2 < U2.
nex());
 
  166   int is = W.
nvol() * i_thread / Nthread;
 
  167   int ns = W.
nvol() * (i_thread + 1) / Nthread - is;
 
  169   double       *g  = W.
ptr(0, is, ex);
 
  170   const double *g1 = U1.
ptr(0, is, ex1);
 
  171   const double *g2 = U2.
ptr(0, is, ex2);
 
  173   const int Nc  = W.
nc();
 
  174   const int Nc2 = 2 * Nc;
 
  175   const int Ndf = Nc2 * Nc;
 
  177   for (
int site = 0; site < ns; ++site) {
 
  180     for (
int ic2 = 0; ic2 < Nc; ++ic2) {
 
  181       for (
int ic1 = 0; ic1 < Nc; ++ic1) {
 
  182         int jg2 = ic2 * 2 + ig;
 
  183         int jg1 = ic1 * 2 + ig;
 
  184         g[2 * ic2 + Nc2 * ic1 + ig]     = mult_Gdn_r(&g1[jg1], &g2[jg2], Nc);
 
  185         g[2 * ic2 + 1 + Nc2 * ic1 + ig] = mult_Gdn_i(&g1[jg1], &g2[jg2], Nc);
 
  197   assert(ex < W.
nex());
 
  198   assert(ex1 < U1.
nex());
 
  199   assert(ex2 < U2.
nex());
 
  208   int is = W.
nvol() * i_thread / Nthread;
 
  209   int ns = W.
nvol() * (i_thread + 1) / Nthread - is;
 
  211   double       *g  = W.
ptr(0, is, ex);
 
  212   const double *g1 = U1.
ptr(0, is, ex1);
 
  213   const double *g2 = U2.
ptr(0, is, ex2);
 
  215   const int Nc  = W.
nc();
 
  216   const int Nc2 = 2 * Nc;
 
  217   const int Ndf = Nc2 * Nc;
 
  219   for (
int site = 0; site < ns; ++site) {
 
  222     for (
int ic2 = 0; ic2 < Nc; ++ic2) {
 
  223       for (
int ic1 = 0; ic1 < Nc; ++ic1) {
 
  224         int jg2 = ic2 * Nc2 + ig;
 
  225         int jg1 = ic1 * Nc2 + ig;
 
  226         g[2 * ic2 + Nc2 * ic1 + ig]     = mult_Gnd_r(&g1[jg1], &g2[jg2], Nc);
 
  227         g[2 * ic2 + 1 + Nc2 * ic1 + ig] = mult_Gnd_i(&g1[jg1], &g2[jg2], Nc);
 
  239   assert(ex < W.
nex());
 
  240   assert(ex1 < U1.
nex());
 
  241   assert(ex2 < U2.
nex());
 
  250   int is = W.
nvol() * i_thread / Nthread;
 
  251   int ns = W.
nvol() * (i_thread + 1) / Nthread - is;
 
  253   double       *g  = W.
ptr(0, is, ex);
 
  254   const double *g1 = U1.
ptr(0, is, ex1);
 
  255   const double *g2 = U2.
ptr(0, is, ex2);
 
  257   const int Nc  = W.
nc();
 
  258   const int Nc2 = 2 * Nc;
 
  259   const int Ndf = Nc2 * Nc;
 
  261   for (
int site = 0; site < ns; ++site) {
 
  264     for (
int ic2 = 0; ic2 < Nc; ++ic2) {
 
  265       for (
int ic1 = 0; ic1 < Nc; ++ic1) {
 
  266         int jg2 = ic2 * Nc2 + ig;
 
  267         int jg1 = ic1 * 2 + ig;
 
  268         g[2 * ic2 + Nc2 * ic1 + ig]     = mult_Gdd_r(&g1[jg1], &g2[jg2], Nc);
 
  269         g[2 * ic2 + 1 + Nc2 * ic1 + ig] = mult_Gdd_i(&g1[jg1], &g2[jg2], Nc);
 
  282   assert(ex < W.
nex());
 
  283   assert(ex1 < U1.
nex());
 
  284   assert(ex2 < U2.
nex());
 
  293   int is = W.
nvol() * i_thread / Nthread;
 
  294   int ns = W.
nvol() * (i_thread + 1) / Nthread - is;
 
  296   double       *g  = W.
ptr(0, is, ex);
 
  297   const double *g1 = U1.
ptr(0, is, ex1);
 
  298   const double *g2 = U2.
ptr(0, is, ex2);
 
  300   const int Nc  = W.
nc();
 
  301   const int Nc2 = 2 * Nc;
 
  302   const int Ndf = Nc2 * Nc;
 
  304   for (
int site = 0; site < ns; ++site) {
 
  307     for (
int ic2 = 0; ic2 < Nc; ++ic2) {
 
  308       for (
int ic1 = 0; ic1 < Nc; ++ic1) {
 
  309         int jg2 = ic2 * 2 + ig;
 
  310         int jg1 = ic1 * Nc2 + ig;
 
  311         g[2 * ic2 + Nc2 * ic1 + ig] +=
 
  312           ff * mult_Gnn_r(&g1[jg1], &g2[jg2], Nc);
 
  313         g[2 * ic2 + 1 + Nc2 * ic1 + ig] +=
 
  314           ff * mult_Gnn_i(&g1[jg1], &g2[jg2], Nc);
 
  327   assert(ex < W.
nex());
 
  328   assert(ex1 < U1.
nex());
 
  329   assert(ex2 < U2.
nex());
 
  338   int is = W.
nvol() * i_thread / Nthread;
 
  339   int ns = W.
nvol() * (i_thread + 1) / Nthread - is;
 
  341   double       *g  = W.
ptr(0, is, ex);
 
  342   const double *g1 = U1.
ptr(0, is, ex1);
 
  343   const double *g2 = U2.
ptr(0, is, ex2);
 
  345   const int Nc  = W.
nc();
 
  346   const int Nc2 = 2 * Nc;
 
  347   const int Ndf = Nc2 * Nc;
 
  349   for (
int site = 0; site < ns; ++site) {
 
  352     for (
int ic2 = 0; ic2 < Nc; ++ic2) {
 
  353       for (
int ic1 = 0; ic1 < Nc; ++ic1) {
 
  354         int jg2 = ic2 * 2 + ig;
 
  355         int jg1 = ic1 * 2 + ig;
 
  356         g[2 * ic2 + Nc2 * ic1 + ig] +=
 
  357           ff * mult_Gdn_r(&g1[jg1], &g2[jg2], Nc);
 
  358         g[2 * ic2 + 1 + Nc2 * ic1 + ig] +=
 
  359           ff * mult_Gdn_i(&g1[jg1], &g2[jg2], Nc);
 
  372   assert(ex < W.
nex());
 
  373   assert(ex1 < U1.
nex());
 
  374   assert(ex2 < U2.
nex());
 
  383   int is = W.
nvol() * i_thread / Nthread;
 
  384   int ns = W.
nvol() * (i_thread + 1) / Nthread - is;
 
  386   double       *g  = W.
ptr(0, is, ex);
 
  387   const double *g1 = U1.
ptr(0, is, ex1);
 
  388   const double *g2 = U2.
ptr(0, is, ex2);
 
  390   const int Nc  = W.
nc();
 
  391   const int Nc2 = 2 * Nc;
 
  392   const int Ndf = Nc2 * Nc;
 
  394   for (
int site = 0; site < ns; ++site) {
 
  397     for (
int ic2 = 0; ic2 < Nc; ++ic2) {
 
  398       for (
int ic1 = 0; ic1 < Nc; ++ic1) {
 
  399         int jg2 = ic2 * Nc2 + ig;
 
  400         int jg1 = ic1 * Nc2 + ig;
 
  401         g[2 * ic2 + Nc2 * ic1 + ig] +=
 
  402           ff * mult_Gnd_r(&g1[jg1], &g2[jg2], Nc);
 
  403         g[2 * ic2 + 1 + Nc2 * ic1 + ig]
 
  404           += ff * mult_Gnd_i(&g1[jg1], &g2[jg2], Nc);
 
  417   assert(ex < W.
nex());
 
  418   assert(ex1 < U1.
nex());
 
  419   assert(ex2 < U2.
nex());
 
  428   int is = W.
nvol() * i_thread / Nthread;
 
  429   int ns = W.
nvol() * (i_thread + 1) / Nthread - is;
 
  431   double       *g  = W.
ptr(0, is, ex);
 
  432   const double *g1 = U1.
ptr(0, is, ex1);
 
  433   const double *g2 = U2.
ptr(0, is, ex2);
 
  435   const int Nc  = W.
nc();
 
  436   const int Nc2 = 2 * Nc;
 
  437   const int Ndf = Nc2 * Nc;
 
  439   for (
int site = 0; site < ns; ++site) {
 
  442     for (
int ic2 = 0; ic2 < Nc; ++ic2) {
 
  443       for (
int ic1 = 0; ic1 < Nc; ++ic1) {
 
  444         int jg2 = ic2 * Nc2 + ig;
 
  445         int jg1 = ic1 * 2 + ig;
 
  446         g[2 * ic2 + Nc2 * ic1 + ig] +=
 
  447           ff * mult_Gdd_r(&g1[jg1], &g2[jg2], Nc);
 
  448         g[2 * ic2 + 1 + Nc2 * ic1 + ig] +=
 
  449           ff * mult_Gdd_i(&g1[jg1], &g2[jg2], Nc);
 
  459   assert(ex < W.
nex());
 
  466   int is = W.
nvol() * i_thread / Nthread;
 
  467   int ns = W.
nvol() * (i_thread + 1) / Nthread - is;
 
  469   double *g = W.
ptr(0, is, ex);
 
  471   const int Nc  = W.
nc();
 
  472   const int Ndf = 2 * Nc * Nc;
 
  474   for (
int site = 0; site < ns; ++site) {
 
  477     for (
int a = 0; a < Nc; ++a) {
 
  478       for (
int b = a + 1; b < Nc; ++b) {
 
  479         double re = g[2 * (Nc * a + b) + ig] - g[2 * (Nc * b + a) + ig];
 
  480         double im = g[2 * (Nc * a + b) + 1 + ig] + g[2 * (Nc * b + a) + 1 + ig];
 
  482         g[2 * (Nc * a + b) + ig]     = 0.5 * re;
 
  483         g[2 * (Nc * a + b) + 1 + ig] = 0.5 * im;
 
  485         g[2 * (Nc * b + a) + ig]     = -0.5 * re;
 
  486         g[2 * (Nc * b + a) + 1 + ig] = 0.5 * im;
 
  490     for (
int cc = 0; cc < Nc; ++cc) {
 
  491       tr += g[2 * (Nc * cc + cc) + 1 + ig];
 
  494     for (
int cc = 0; cc < Nc; ++cc) {
 
  495       g[2 * (Nc * cc + cc) + ig]      = 0.0;
 
  496       g[2 * (Nc * cc + cc) + 1 + ig] -= tr;
 
  505   assert(ex < W.
nex());
 
  512   int is = W.
nvol() * i_thread / Nthread;
 
  513   int ns = W.
nvol() * (i_thread + 1) / Nthread - is;
 
  515   double *g = W.
ptr(0, is, ex);
 
  517   const int Nc  = W.
nc();
 
  518   const int Ndf = 2 * Nc * Nc;
 
  520   for (
int site = 0; site < ns; ++site) {
 
  523     for (
int a = 0; a < Nc; ++a) {
 
  524       for (
int b = a; b < Nc; ++b) {
 
  525         double re = g[2 * (Nc * a + b) + ig] - g[2 * (Nc * b + a) + ig];
 
  526         double im = g[2 * (Nc * a + b) + 1 + ig] + g[2 * (Nc * b + a) + 1 + ig];
 
  528         g[2 * (Nc * a + b) + ig]     = 0.5 * re;
 
  529         g[2 * (Nc * a + b) + 1 + ig] = 0.5 * im;
 
  531         g[2 * (Nc * b + a) + ig]     = -0.5 * re;
 
  532         g[2 * (Nc * b + a) + 1 + ig] = 0.5 * im;
 
  552   for (
int ex = 0; ex < Nex; ++ex) {
 
  553     for (
int site = 0; site < Nvol; ++site) {
 
  554       u0 = U.
mat(site, ex);
 
  555       v0 = iP.
mat(site, ex);
 
const double * ptr(const int jin, const int site, const int jex) const 
 
void at_Field_G(Field_G &W, const int ex)
 
void check()
check of assumptions for performance implementation. 
 
void ah_Field_G(Field_G &W, const int ex)
 
virtual void gauss_lex_global(Field &)
gaussian random number defined on global lattice. 
 
void set_random(RandomNumbers *rand)
 
void multadd_Field_Gdd(Field_G &W, const int ex, const Field_G &U1, int ex1, const Field_G &U2, int ex2, const double ff)
 
void multadd_Field_Gnn(Field_G &W, const int ex, const Field_G &U1, int ex1, const Field_G &U2, int ex2, const double ff)
 
void multadd_Field_Gdn(Field_G &W, const int ex, const Field_G &U1, int ex1, const Field_G &U2, int ex2, const double ff)
 
void multadd_Field_Gnd(Field_G &W, const int ex, const Field_G &U1, int ex1, const Field_G &U2, int ex2, const double ff)
 
void mult_Field_Gnn(Field_G &W, const int ex, const Field_G &U1, int ex1, const Field_G &U2, int ex2)
 
void mult_Field_Gnd(Field_G &W, const int ex, const Field_G &U1, int ex1, const Field_G &U2, int ex2)
 
Base class of random number generators. 
 
void mult_Field_Gdd(Field_G &W, const int ex, const Field_G &U1, int ex1, const Field_G &U2, int ex2)
 
Mat_SU_N mat_exp(const double alpha, const Mat_SU_N &iv, const Mat_SU_N &u, const int Nprec)
 
void set_mat(const int site, const int mn, const Mat_SU_N &U)
 
void mult_exp_Field_G(Field_G &W, const double alpha, const Field_G &iP, const Field_G &U, const int Nprec)
 
Mat_SU_N mat(const int site, const int mn=0) const 
 
void mult_Field_Gdn(Field_G &W, const int ex, const Field_G &U1, int ex1, const Field_G &U2, int ex2)