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);
141 for (
int mu = 0; mu <
m_Ndim; ++mu) {
165 params_solver.
set_string(
"solver_type",
"CG");
166 params_solver.
set_int(
"maximum_number_of_iteration", 1000);
167 params_solver.
set_int(
"maximum_number_of_restart", 40);
168 params_solver.
set_double(
"convergence_criterion_squared", 1.0e-30);
170 params_solver.
set_string(
"verbose_level",
"Crucial");
183 const int Niter = 100;
184 const int Nrestart = 40;
185 const double Stopping_condition = 1.0e-30;
194 for (
int ispin = 0; ispin <
m_Nd; ++ispin) {
195 for (
int icolor = 0; icolor <
m_Nc; ++icolor) {
196 int spin_color = icolor + m_Nc * ispin;
198 for (
int isite = 0; isite <
m_Nvol2; ++isite) {
199 w.
set_ri(icolor, ispin, isite, 0, 1, 0);
207 solver->
solve(w2, w, Nconv, diff);
212 solver->
solve(w2, w, Nconv, diff);
221 for (
int ics = 0; ics <
m_Nc *
m_Nd; ++ics) {
222 for (
int site = 0; site <
m_Nvol2; ++site) {
223 for (
int id = 0;
id <
m_Nd; ++id) {
224 for (
int ic = 0; ic <
m_Nc; ++ic) {
255 const Field& f,
const int ieo)
259 }
else if (
m_repr ==
"Chiral") {
267 const Field& f,
const int ieo)
271 const double *
v1 = f.
ptr(0);
272 double *
v2 = v.
ptr(0);
277 }
else if (ieo == 1) {
294 int ns =
m_Nvol2 * (ith + 1) / nth;
297 for (
int site = is; site < ns; ++site) {
298 for (
int icd = 0; icd <
m_Nc * Nd2; ++icd) {
299 int iv2 = 2 * icd +
m_NinF * site;
302 for (
int jd = 0; jd <
m_Nd; ++jd) {
304 int iv = jcd +
m_NinF * site;
306 v2[iv2] += mult_uv_r(&csw_inv[ig], &v1[iv],
m_Nc);
307 v2[iv2 + 1] += mult_uv_i(&csw_inv[ig], &v1[iv],
m_Nc);
311 for (
int icd = 0; icd <
m_Nc * Nd2; ++icd) {
312 int iv2 = 2 * (icd +
m_Nc * Nd2) +
m_NinF * site;
315 for (
int jd = 0; jd <
m_Nd; ++jd) {
316 int jd2 = (jd + Nd2) % m_Nd;
317 int iv = Nvc * jd +
m_NinF * site;
319 v2[iv2] += mult_uv_r(&csw_inv[ig], &v1[iv],
m_Nc);
320 v2[iv2 + 1] += mult_uv_i(&csw_inv[ig], &v1[iv],
m_Nc);
330 const Field& f,
const int ieo)
334 const double *
v1 = f.
ptr(0);
335 double *
v2 = v.
ptr(0);
340 }
else if (ieo == 1) {
357 int ns =
m_Nvol2 * (ith + 1) / nth;
359 for (
int site = is; site < ns; ++site) {
360 for (
int icd = 0; icd <
m_Nc *
m_Nd / 2; ++icd) {
361 int iv2 = 2 * icd +
m_NinF * site;
365 for (
int jd = 0; jd <
m_Nd / 2; ++jd) {
367 int iv = jcd +
m_NinF * site;
369 v2[iv2] += mult_uv_r(&csw_inv[ig], &v1[iv],
m_Nc);
370 v2[iv2 + 1] += mult_uv_i(&csw_inv[ig], &v1[iv],
m_Nc);
375 int iv2 = 2 * icd +
m_NinF * site;
379 for (
int jd = m_Nd / 2; jd <
m_Nd; ++jd) {
381 int iv = jcd +
m_NinF * site;
383 v2[iv2] += mult_uv_r(&csw_inv[ig], &v1[iv],
m_Nc);
384 v2[iv2 + 1] += mult_uv_i(&csw_inv[ig], &v1[iv],
m_Nc);
397 for (
int ispin = 0; ispin < m_Nd / 2; ++ispin) {
398 for (
int icolor = 0; icolor <
m_Nc; ++icolor) {
399 int ics = icolor + ispin *
m_Nc;
400 for (
int jspin = 0; jspin <
m_Nd; ++jspin) {
401 int js2 = (jspin + m_Nd / 2) % m_Nd;
402 for (
int jcolor = 0; jcolor <
m_Nc; ++jcolor) {
403 int cs1 = jcolor + m_Nc * (jspin + m_Nd * ics);
404 int cs2 = jcolor + m_Nc * (jspin + m_Nd * (ics + m_Nc * m_Nd / 2));
405 int cc = jcolor + icolor *
m_Nc;
406 int ss1 = jspin + ispin *
m_Nd;
407 int ss2 = js2 + ispin *
m_Nd;
409 matrix[2 * cs1] =
m_T.
cmp_r(cc, site, ss1);
410 matrix[2 * cs1 + 1] =
m_T.
cmp_i(cc, site, ss1);
412 matrix[2 * cs2] =
m_T.
cmp_r(cc, site, ss2);
413 matrix[2 * cs2 + 1] =
m_T.
cmp_i(cc, site, ss2);
428 }
else if (
m_repr ==
"Chiral") {
442 const double *fp = f.
ptr(0);
443 double *vp = v.
ptr(0);
449 int ns =
m_Nvol2 * (ith + 1) / nth;
456 for (
int site = is; site < ns; ++site) {
457 for (
int id = 0;
id < Nd2; ++id) {
458 for (
int ic = 0; ic <
m_Nc; ++ic) {
459 int icd = ic + m_Nc * id;
461 int iv2 = 2 * icd + NinF * site;
464 for (
int jd = 0; jd <
m_Nd; ++jd) {
465 int iv = Nvc * jd + NinF * site;
466 int ig = Nvc * ic + NinG * (site +
m_Nvol * (
id * m_Nd + jd));
467 vp[iv2] += mult_uv_r(&tp[ig], &fp[iv], m_Nc);
468 vp[iv2 + 1] += mult_uv_i(&tp[ig], &fp[iv], m_Nc);
474 for (
int jd = 0; jd <
m_Nd; ++jd) {
475 int jd2 = (2 + jd) % m_Nd;
476 int iv = Nvc * jd + NinF * site;
477 int ig = Nvc * ic + NinG * (site +
m_Nvol * (
id * m_Nd + jd2));
478 vp[iv2] += mult_uv_r(&tp[ig], &fp[iv], m_Nc);
479 vp[iv2 + 1] += mult_uv_i(&tp[ig], &fp[iv], m_Nc);
491 const double *fp = f.
ptr(0);
492 double *vp = v.
ptr(0);
498 int ns =
m_Nvol2 * (ith + 1) / nth;
505 for (
int site = is; site < ns; ++site) {
506 for (
int id = 0;
id < Nd2; ++id) {
507 for (
int ic = 0; ic <
m_Nc; ++ic) {
508 int icd = ic + m_Nc * id;
510 int iv2 = 2 * icd + NinF * site;
513 for (
int jd = 0; jd < Nd2; ++jd) {
514 int iv = Nvc * jd + NinF * site;
515 int ig = Nvc * ic + NinG * (site +
m_Nvol * (
id * Nd2 + jd));
516 vp[iv2] += mult_uv_r(&tp[ig], &fp[iv], m_Nc);
517 vp[iv2 + 1] += mult_uv_i(&tp[ig], &fp[iv], m_Nc);
523 for (
int jd = 0; jd < Nd2; ++jd) {
524 int iv = Nvc * (Nd2 + jd) + NinF * site;
525 int ig = Nvc * ic + NinG * (site +
m_Nvol * (m_Nd +
id * Nd2 + jd));
526 vp[iv2] += mult_uv_r(&tp[ig], &fp[iv], m_Nc);
527 vp[iv2 + 1] += mult_uv_i(&tp[ig], &fp[iv], m_Nc);
538 const int mu,
const int nu)
550 }
else if (
m_repr ==
"Chiral") {
553 vout.
crucial(
m_vl,
"Error at %s: unsupported gamma matrix repr. %s.\n",
708 const int mu,
const int nu)
721 for (
int site = 0; site <
m_Nvol; ++site) {
722 w.set_mat(site, 0, Umu.
mat(site) * Cup.mat_dag(site));
725 for (
int site = 0; site <
m_Nvol; ++site) {
731 for (
int site = 0; site <
m_Nvol; ++site) {
732 v.set_mat(site, 0, Cup.mat_dag(site) * Umu.
mat(site));
735 for (
int site = 0; site <
m_Nvol; ++site) {
745 for (
int site = 0; site <
m_Nvol; ++site) {
746 Fst.
set_mat(site, 0, w.mat(site).ah());
761 assert(tr_sigma_inv.
nex() == 1);
771 for (
int isite = 0; isite <
m_Nvol; ++isite) {
772 for (
int ispin = 0; ispin <
m_Nd; ++ispin) {
773 for (
int icolor = 0; icolor <
m_Nc; ++icolor) {
774 v = sigma_inv.
vec(ispin, isite, icolor + m_Nc * ispin);
775 for (
int jcolor = 0; jcolor <
m_Nc; ++jcolor) {
776 int cc = icolor + m_Nc * jcolor;
777 tr_sigma_inv.
set_r(cc, isite, 0, v.
r(jcolor));
778 tr_sigma_inv.
set_i(cc, isite, 0, v.
i(jcolor));
794 double flop_site = 0.0;
798 }
else if (
m_repr ==
"Chiral") {
802 double flop = flop_site *
static_cast<double>(Lvol / 2);
void scal(Field &x, const double a)
scal(x, a): x = a * x
static const std::string class_name
void set_config(Field *Ueo)
setting pointer to the gauge configuration.
void init(std::string repr)
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.
void mult_csw_inv_chiral(Field &, const Field &, const int ieo)
const double * ptr(const int jin, const int site, const int jex) const
void mult_csw_inv(Field &, const Field &, const int ieo)
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.
void general(const char *format,...)
GammaMatrix get_GM(GMspecies spec)
void set_int(const string &key, const int value)
std::vector< int > m_boundary
void D_dirac(Field &v, const Field &f, const int ieo)
explicit implementation for Dirac representation (for Imp-version).
Container of Field-type object.
int fetch_double(const string &key, double &value) const
double cmp_i(const int cc, const int site, const int mn=0) const
void set_csw_chiral()
explicit implementation for Chiral representation (for Imp-version).
void mult_isigma(Field_F &, const Field_F &, const int mu, const int nu)
Standard Conjugate Gradient solver algorithm.
std::vector< GammaMatrix > m_SG
void mult_csw_inv_dirac(Field &, const Field &, const int ieo)
static int get_thread_id()
returns thread id.
Wilson-type fermion field.
virtual void set_parameters(const Parameters ¶ms)=0
void set_csw_dirac()
explicit implementation for Dirac representation (for Imp-version).
void set_string(const string &key, const string &value)
Field_G m_T
m_T = 1 - kappa c_SW sigma F / 2
static double epsilon_criterion2()
void reset(int Nvol, int Nex)
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_ri(const int cc, const int s, const int site, const int e, const double re, const double im)
Bridge::VerboseLevel m_vl
void set_i(const int cc, const int site, const int mn, const double im)
std::vector< double > csmatrix(const int &)
void set_r(const int cc, const int site, const int mn, const double re)
void D_chiral(Field &v, const Field &f, const int ieo)
explicit implementation for Chiral representation (for Imp-version).
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
void crucial(const char *format,...)
void set_double(const string &key, const double value)
int sg_index(int mu, int nu)
std::vector< GammaMatrix > m_GM
Gamma Matrix and Sigma_{mu,nu} = -i [Gamma_mu, Gamma_nu] /2.
void reverseField(Field &lex, const Field &eo)
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_mode(std::string mode)
setting the mode of multiplication if necessary. Default implementation here is just to avoid irrelev...
double cmp_r(const int cc, const int site, const int mn=0) const
void set_fieldstrength(Field_G &, const int, const int)
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)
void set_parameters(const Parameters ¶ms)
Mat_SU_N mat(const int site, const int mn=0) const
double flop_count()
retuns number of floating point number operations.
virtual void solve(Field &solution, const Field &source, int &Nconv, double &diff)=0
static VerboseLevel set_verbose_level(const std::string &str)
void trSigmaInv(Field_G &, const int mu, const int nu)
void D(Field &v, const Field &f, const int ieo)
double cmp_r(const int cc, const int s, const int site, const int e=0) const