Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
action_F_Standard_SF.cpp
Go to the documentation of this file.
1 
14 #include "action_F_Standard_SF.h"
15 
16 //====================================================================
17 void Action_F_Standard_SF::set_parameters(const Parameters& params)
18 {
19  Bridge::VerboseLevel vl = params.get_VerboseLevel();
20 
21  set_parameter_verboselevel(vl);
22 }
23 
24 
25 //====================================================================
26 void Action_F_Standard_SF::set_parameters()
27 {
28  int Nc = CommonParameters::Nc();
29  int Nvol = CommonParameters::Nvol();
30  int Ndim = CommonParameters::Ndim();
31  int NinG = 2 * Nc * Nc;
32 
33  vout.general(m_vl, "Action_F_Standard_SF:\n");
34 
35  m_force.reset(NinG, Nvol, Ndim);
36 
37  int Niter = 20000;
38  double Stop_cond = 1.0e-24;
39 
40 
41  string str_solver_type = "CG";
42  m_solver = Solver::New(str_solver_type, m_fopr);
43  m_solver->set_parameters(Niter, Stop_cond);
44 
45  // link variable update flag
46  m_status_linkv = 0;
47 }
48 
49 
50 //====================================================================
51 double Action_F_Standard_SF::langevin(RandomNumbers *rand)
52 {
53  int Nvol = CommonParameters::Nvol();
54  int Ndim = CommonParameters::Ndim();
55 
56  int NinF = m_fopr->field_nin();
57  int NvolF = m_fopr->field_nvol();
58  int NexF = m_fopr->field_nex();
59  int size_psf = NinF * NvolF * NexF * CommonParameters::NPE();
60 
61  assert(NvolF == Nvol);
62  m_psf.reset(NinF, NvolF, NexF);
63 
64  vout.general(m_vl, " Action_F_Standard_SF: %s\n", m_label.c_str());
65 
66  Field xi(NinF, NvolF, NexF);
67  rand->gauss_lex_global(xi);
68 
69  m_fopr->set_config(m_U);
70  m_fopr->set_mode("Ddag");
71  m_psf = m_fopr->mult(xi);
72 
73  // set_boundary_zero(xi);
74  Field_F_SF setzero;
75  setzero.set_boundary_zero(xi);
76  double xi2 = xi.norm();
77  double H_psf = xi2 * xi2;
78 
79  vout.general(m_vl, " H_Fstandard = %18.8f\n", H_psf);
80  vout.general(m_vl, " H_F/dof = %18.8f\n", H_psf / size_psf);
81 
82  return H_psf;
83 }
84 
85 
86 //====================================================================
87 double Action_F_Standard_SF::calcH()
88 {
89  int Nvol = CommonParameters::Nvol();
90  int Ndim = CommonParameters::Ndim();
91 
92  int NinF = m_fopr->field_nin();
93  int NvolF = m_fopr->field_nvol();
94  int NexF = m_fopr->field_nex();
95  int size_psf = NinF * NvolF * NexF * CommonParameters::NPE();
96 
97  Field v1(NinF, NvolF, NexF);
98 
99  vout.general(m_vl, " Action_F_Standard_SF: %s\n", m_label.c_str());
100 
101  m_fopr->set_config(m_U);
102  m_fopr->set_mode("DdagD");
103 
104  int Nconv;
105  double diff;
106  m_solver->solve(v1, m_psf, Nconv, diff);
107 
108  vout.general(m_vl, " Nconv = %d diff = %.8e\n", Nconv, diff);
109 
110  Field_F_SF setzero;
111  setzero.set_boundary_zero(v1);
112  setzero.set_boundary_zero(m_psf);
113 
114  double H_psf = v1 * m_psf;
115 
116  vout.general(m_vl, " H_Fstandard = %18.8f\n", H_psf);
117  vout.general(m_vl, " H_F/dof = %18.8f\n", H_psf / size_psf);
118 
119  return H_psf;
120 }
121 
122 
123 //====================================================================
124 const Field Action_F_Standard_SF::force()
125 {
126  if (m_status_linkv == 0) {
127  int Nin = m_U->nin();
128  int Nvol = m_U->nvol();
129  int Nex = m_U->nex();
130  int Nc = CommonParameters::Nc();
131  int Ndim = CommonParameters::Ndim();
132 
133  int NinF = m_fopr->field_nin();
134  int NvolF = m_fopr->field_nvol();
135  int NexF = m_fopr->field_nex();
136  Field eta(NinF, NvolF, NexF);
137 
138  vout.general(m_vl, " Action_F_Standard_SF: %s\n", m_label.c_str());
139 
140  //- fermion inversion for smeared gauge field
141  m_fopr->set_config(m_U);
142  m_fopr->set_mode("DdagD");
143 
144  int Nconv;
145  double diff;
146  m_solver->solve(eta, m_psf, Nconv, diff);
147 
148  vout.general(m_vl, " Solver: Nconv = %6d diff = %12.6e\n", Nconv, diff);
149 
150  //- force of smeared fermion operator
151  m_fopr_force->set_config(m_U);
152 
153  Field force(Nin, Nvol, Nex);
154  force = m_fopr_force->force_core(eta);
155 
156  Field_G_SF Fboundary(force);
157  Fboundary.set_boundary_spatial_link_zero();
158  m_force = (Field)Fboundary;
159 
160  ++m_status_linkv;
161 
162  double Fave, Fmax, Fdev;
163  m_force.stat(Fave, Fmax, Fdev);
164  vout.general(m_vl, " Fstandard_ave = %12.6f Fstandard_max = %12.6f Fstandard_dev = %12.6f\n",
165  Fave, Fmax, Fdev);
166 
167  return m_force;
168  } else {
169  vout.general(m_vl, " Fstandard returns previous force.\n");
170  return m_force;
171  }
172 }
173 
174 
175 //====================================================================
176 
181 /*
182 void Action_F_Standard_SF::set_boundary_zero(Field& f){
183  if(comm->ipe(3)==0){
184  for(int site = 0; site < Svol; ++site){
185  for(int s = 0; s < m_Nd; ++s){
186  for(int cc = 0; cc < m_Nc2; ++cc){
187  f.set(cc+m_Nc2*s, site, 0, 0.0);
188  }
189  }
190  }
191  }
192 }
193 */
194 
195 //====================================================================
196 //============================================================END=====