21 #if defined USE_GROUP_SU3
22 #include "fopr_Wilson_impl_SU3.inc"
23 #elif defined USE_GROUP_SU2
24 #include "fopr_Wilson_impl_SU2.inc"
25 #elif defined USE_GROUP_SU_N
26 #include "fopr_Wilson_impl_SU_N.inc"
99 const string str_vlevel = params.
get_string(
"verbose_level");
123 const std::vector<int> bc)
129 for (
int mu = 0; mu <
m_Ndim; ++mu) {
135 assert(bc.size() ==
m_Ndim);
140 assert(bc.size() ==
m_Ndim);
142 for (
int mu = 0; mu <
m_Ndim; ++mu) {
166 params_solver.
set_string(
"solver_type",
"CG");
167 params_solver.
set_int(
"maximum_number_of_iteration", 100);
168 params_solver.
set_int(
"maximum_number_of_restart", 40);
169 params_solver.
set_double(
"convergence_criterion_squared", 1.0e-30);
171 params_solver.
set_string(
"verbose_level",
"Crucial");
184 const int Niter = 100;
185 const int Nrestart = 40;
186 const double Stopping_condition = 1.0e-30;
195 for (
int ispin = 0; ispin <
m_Nd; ++ispin) {
196 for (
int icolor = 0; icolor <
m_Nc; ++icolor) {
197 int spin_color = icolor + m_Nc * ispin;
199 for (
int isite = 0; isite <
m_Nvol2; ++isite) {
200 w.
set_ri(icolor, ispin, isite, 0, 1, 0);
208 solver->
solve(w2, w, Nconv, diff);
213 solver->
solve(w2, w, Nconv, diff);
222 for (
int ics = 0; ics <
m_Nc *
m_Nd; ++ics) {
223 for (
int site = 0; site <
m_Nvol2; ++site) {
224 for (
int id = 0;
id <
m_Nd; ++id) {
225 for (
int ic = 0; ic <
m_Nc; ++ic) {
242 const Field& f,
const int ieo)
246 }
else if (
m_repr ==
"Chiral") {
254 const Field& f,
const int ieo)
258 const double *
v1 = f.
ptr(0);
259 double *
v2 = v.
ptr(0);
265 }
else if (ieo == 1) {
279 int is =
m_Nvol2 * i_thread / Nthread;
280 int ns =
m_Nvol2 * (i_thread + 1) / Nthread;
283 for (
int site = is; site < ns; ++site) {
284 for (
int icd = 0; icd <
m_Nc * Nd2; ++icd) {
285 int iv2 = 2 * icd +
m_NinF * site;
288 for (
int jd = 0; jd <
m_Nd; ++jd) {
290 int iv = jcd +
m_NinF * site;
292 v2[iv2] += mult_uv_r(&csw_inv[ig], &v1[iv],
m_Nc);
293 v2[iv2 + 1] += mult_uv_i(&csw_inv[ig], &v1[iv],
m_Nc);
297 for (
int icd = 0; icd <
m_Nc * Nd2; ++icd) {
298 int iv2 = 2 * (icd +
m_Nc * Nd2) +
m_NinF * site;
301 for (
int jd = 0; jd <
m_Nd; ++jd) {
302 int jd2 = (jd + Nd2) % m_Nd;
303 int iv = Nvc * jd +
m_NinF * site;
305 v2[iv2] += mult_uv_r(&csw_inv[ig], &v1[iv],
m_Nc);
306 v2[iv2 + 1] += mult_uv_i(&csw_inv[ig], &v1[iv],
m_Nc);
316 const Field& f,
const int ieo)
320 const double *
v1 = f.
ptr(0);
321 double *
v2 = v.
ptr(0);
327 }
else if (ieo == 1) {
341 int is =
m_Nvol2 * i_thread / Nthread;
342 int ns =
m_Nvol2 * (i_thread + 1) / Nthread;
344 for (
int site = is; site < ns; ++site) {
345 for (
int icd = 0; icd <
m_Nc *
m_Nd / 2; ++icd) {
346 int iv2 = 2 * icd +
m_NinF * site;
350 for (
int jd = 0; jd <
m_Nd / 2; ++jd) {
352 int iv = jcd +
m_NinF * site;
354 v2[iv2] += mult_uv_r(&csw_inv[ig], &v1[iv],
m_Nc);
355 v2[iv2 + 1] += mult_uv_i(&csw_inv[ig], &v1[iv],
m_Nc);
360 int iv2 = 2 * icd +
m_NinF * site;
364 for (
int jd = m_Nd / 2; jd <
m_Nd; ++jd) {
366 int iv = jcd +
m_NinF * site;
368 v2[iv2] += mult_uv_r(&csw_inv[ig], &v1[iv],
m_Nc);
369 v2[iv2 + 1] += mult_uv_i(&csw_inv[ig], &v1[iv],
m_Nc);
382 for (
int ispin = 0; ispin < m_Nd / 2; ++ispin) {
383 for (
int icolor = 0; icolor <
m_Nc; ++icolor) {
384 int ics = icolor + ispin *
m_Nc;
385 for (
int jspin = 0; jspin <
m_Nd; ++jspin) {
386 int js2 = (jspin + m_Nd / 2) % m_Nd;
387 for (
int jcolor = 0; jcolor <
m_Nc; ++jcolor) {
388 int cs1 = jcolor + m_Nc * (jspin + m_Nd * ics);
389 int cs2 = jcolor + m_Nc * (jspin + m_Nd * (ics + m_Nc * m_Nd / 2));
390 int cc = jcolor + icolor *
m_Nc;
391 int ss1 = jspin + ispin *
m_Nd;
392 int ss2 = js2 + ispin *
m_Nd;
394 matrix[2 * cs1] =
m_T.
cmp_r(cc, site, ss1);
395 matrix[2 * cs1 + 1] =
m_T.
cmp_i(cc, site, ss1);
397 matrix[2 * cs2] =
m_T.
cmp_r(cc, site, ss2);
398 matrix[2 * cs2 + 1] =
m_T.
cmp_i(cc, site, ss2);
413 }
else if (
m_repr ==
"Chiral") {
427 const double *fp = f.
ptr(0);
428 double *vp = v.
ptr(0);
433 int is =
m_Nvol2 * i_thread / Nthread;
434 int ns =
m_Nvol2 * (i_thread + 1) / Nthread;
441 for (
int site = is; site < ns; ++site) {
442 for (
int id = 0;
id < Nd2; ++id) {
443 for (
int ic = 0; ic <
m_Nc; ++ic) {
444 int icd = ic + m_Nc * id;
446 int iv2 = 2 * icd + NinF * site;
449 for (
int jd = 0; jd <
m_Nd; ++jd) {
450 int iv = Nvc * jd + NinF * site;
451 int ig = Nvc * ic + NinG * (site +
m_Nvol * (
id * m_Nd + jd));
452 vp[iv2] += mult_uv_r(&tp[ig], &fp[iv], m_Nc);
453 vp[iv2 + 1] += mult_uv_i(&tp[ig], &fp[iv], m_Nc);
459 for (
int jd = 0; jd <
m_Nd; ++jd) {
460 int jd2 = (2 + jd) % m_Nd;
461 int iv = Nvc * jd + NinF * site;
462 int ig = Nvc * ic + NinG * (site +
m_Nvol * (
id * m_Nd + jd2));
463 vp[iv2] += mult_uv_r(&tp[ig], &fp[iv], m_Nc);
464 vp[iv2 + 1] += mult_uv_i(&tp[ig], &fp[iv], m_Nc);
476 const double *fp = f.
ptr(0);
477 double *vp = v.
ptr(0);
482 int is =
m_Nvol2 * i_thread / Nthread;
483 int ns =
m_Nvol2 * (i_thread + 1) / Nthread;
490 for (
int site = is; site < ns; ++site) {
491 for (
int id = 0;
id < Nd2; ++id) {
492 for (
int ic = 0; ic <
m_Nc; ++ic) {
493 int icd = ic + m_Nc * id;
495 int iv2 = 2 * icd + NinF * site;
498 for (
int jd = 0; jd < Nd2; ++jd) {
499 int iv = Nvc * jd + NinF * site;
500 int ig = Nvc * ic + NinG * (site +
m_Nvol * (
id * Nd2 + jd));
501 vp[iv2] += mult_uv_r(&tp[ig], &fp[iv], m_Nc);
502 vp[iv2 + 1] += mult_uv_i(&tp[ig], &fp[iv], m_Nc);
508 for (
int jd = 0; jd < Nd2; ++jd) {
509 int iv = Nvc * (Nd2 + jd) + NinF * site;
510 int ig = Nvc * ic + NinG * (site +
m_Nvol * (m_Nd +
id * Nd2 + jd));
511 vp[iv2] += mult_uv_r(&tp[ig], &fp[iv], m_Nc);
512 vp[iv2 + 1] += mult_uv_i(&tp[ig], &fp[iv], m_Nc);
523 const int mu,
const int nu)
535 }
else if (
m_repr ==
"Chiral") {
538 vout.
crucial(
m_vl,
"Error at %s: unsupported gamma matrix repr. %s.\n",
693 const int mu,
const int nu)
706 for (
int site = 0; site <
m_Nvol; ++site) {
707 w.set_mat(site, 0, Umu.
mat(site) * Cup.mat_dag(site));
710 for (
int site = 0; site <
m_Nvol; ++site) {
716 for (
int site = 0; site <
m_Nvol; ++site) {
717 v.set_mat(site, 0, Cup.mat_dag(site) * Umu.
mat(site));
720 for (
int site = 0; site <
m_Nvol; ++site) {
730 for (
int site = 0; site <
m_Nvol; ++site) {
731 Fst.
set_mat(site, 0, w.mat(site).ah());
746 assert(tr_sigma_inv.
nex() == 1);
755 for (
int isite = 0; isite <
m_Nvol; ++isite) {
756 for (
int ispin = 0; ispin <
m_Nd; ++ispin) {
757 for (
int icolor = 0; icolor <
m_Nc; ++icolor) {
758 v = sigma_inv.
vec(ispin, isite, icolor + m_Nc * ispin);
759 for (
int jcolor = 0; jcolor <
m_Nc; ++jcolor) {
760 int cc = icolor + m_Nc * jcolor;
761 tr_sigma_inv.
set_r(cc, isite, 0, v.
r(jcolor));
762 tr_sigma_inv.
set_i(cc, isite, 0, v.
i(jcolor));
778 double flop_site = 0.0;
782 }
else if (
m_repr ==
"Chiral") {
786 double flop = flop_site *
static_cast<double>(Lvol / 2);
void scal(Field &x, const double a)
scal(x, a): x = a * x
void mult_csw_inv_dirac(Field &, const Field &, const int ieo)
void D(Field &v, const Field &f, const int ieo)
double cmp_i(const int cc, const int s, const int site, const int e=0) const
void detailed(const char *format,...)
static int get_num_threads()
returns available number of threads.
const double * ptr(const int jin, const int site, const int jex) const
virtual void upper(Field_G &, const Field_G &, const int mu, const int nu)=0
constructs upper staple in mu-nu plane.
double r(const int c) const
void set(const int jin, const int site, const int jex, double v)
virtual void lower(Field_G &, const Field_G &, const int mu, const int nu)=0
constructs lower staple in mu-nu plane.
std::vector< GammaMatrix > m_SG
void general(const char *format,...)
GammaMatrix get_GM(GMspecies spec)
void set_int(const string &key, const int value)
void solve(Field &solution, const Field &source, int &Nconv, double &diff)
std::vector< GammaMatrix > m_GM
Gamma Matrix and Sigma_{mu,nu} = -i [Gamma_mu, Gamma_nu] /2.
Container of Field-type object.
int fetch_double(const string &key, double &value) const
void set_parameters(const Parameters ¶ms)
double cmp_i(const int cc, const int site, const int mn=0) const
Standard Conjugate Gradient solver algorithm.
void init(std::string repr)
void set_csw_chiral()
explicit implementation for Chiral representation (for Imp-version).
static int get_thread_id()
returns thread id.
Wilson-type fermion field.
std::vector< int > m_boundary
void trSigmaInv(Field_G &, const int mu, const int nu)
void set_string(const string &key, const string &value)
void set_fieldstrength(Field_G &, const int, const int)
static double epsilon_criterion2()
void reset(int Nvol, int Nex)
void mult_csw_inv_chiral(Field &, const Field &, const int ieo)
void mult_iGM(Field_F &y, const GammaMatrix &gm, const Field_F &x)
gamma matrix multiplication (i is multiplied)
double i(const int c) const
void set_parameters(const Parameters ¶ms)
void set_ri(const int cc, const int s, const int site, const int e, const double re, const double im)
Bridge::VerboseLevel m_vl
double flop_count()
retuns number of floating point number operations.
void set_mode(std::string mode)
setting the mode of multiplication if necessary. Default implementation here is just to avoid irrelev...
void set_i(const int cc, const int site, const int mn, const double im)
std::vector< double > csmatrix(const int &)
void mult_isigma(Field_F &, const Field_F &, const int mu, const int nu)
void set_r(const int cc, const int site, const int mn, const double re)
void mult_csw_inv(Field &, const Field &, const int ieo)
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
void crucial(const char *format,...)
static const std::string class_name
void set_double(const string &key, const double value)
int sg_index(int mu, int nu)
void reverseField(Field &lex, const Field &eo)
void set_config(Field *Ueo)
setting pointer to the gauge configuration.
void forward(Field &, const Field &, const int mu)
Vec_SU_N vec(const int s, const int site, const int e=0) const
Mat_SU_N mat_dag(const int site, const int mn=0) const
void set_csw_dirac()
explicit implementation for Dirac representation (for Imp-version).
double cmp_r(const int cc, const int site, const int mn=0) const
Field_G m_T
m_T = 1 - kappa c_SW sigma F / 2
void D_dirac(Field &v, const Field &f, const int ieo)
explicit implementation for Dirac representation (for Imp-version).
void setpart_ex(int ex, const Field &w, int exw)
void set_parameter_verboselevel(const Bridge::VerboseLevel vl)
string get_string(const string &key) const
int fetch_int_vector(const string &key, vector< int > &value) const
void set_mat(const int site, const int mn, const Mat_SU_N &U)
Mat_SU_N mat(const int site, const int mn=0) const
static VerboseLevel set_verbose_level(const std::string &str)
void D_chiral(Field &v, const Field &f, const int ieo)
explicit implementation for Chiral representation (for Imp-version).
double cmp_r(const int cc, const int s, const int site, const int e=0) const