Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
action_F_Ratio.cpp
Go to the documentation of this file.
1 
14 #include "action_F_Ratio.h"
15 
16 //====================================================================
17 void Action_F_Ratio::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_Ratio::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_Ratio:\n");
34 
35  m_force.reset(NinG, Nvol, Ndim);
36 
37  int Niter = 2000;
38  double Stop_cond_MD = 1.0e-24;
39  double Stop_cond_H = 1.0e-24;
40 
41  string str_solver_type = "CG";
42  m_solver_H_prec = Solver::New(str_solver_type, m_fopr_prec);
43  m_solver_MD = Solver::New(str_solver_type, m_fopr);
44  m_solver_H = Solver::New(str_solver_type, m_fopr);
45 
46  m_solver_H_prec->set_parameters(Niter, Stop_cond_H);
47  m_solver_MD->set_parameters(Niter, Stop_cond_MD);
48  m_solver_H->set_parameters(Niter, Stop_cond_H);
49 
50  //- link variable update flag
51  m_status_linkv = 0;
52 }
53 
54 
55 //====================================================================
56 double Action_F_Ratio::langevin(RandomNumbers *rand)
57 {
58  int Nvol = CommonParameters::Nvol();
59  int Ndim = CommonParameters::Ndim();
60 
61  int NinF = m_fopr_prec->field_nin();
62  int NvolF = m_fopr_prec->field_nvol();
63  int NexF = m_fopr_prec->field_nex();
64  int size_psf = NinF * NvolF * NexF * CommonParameters::NPE();
65 
66  assert(NvolF == Nvol);
67  m_psf.reset(NinF, NvolF, NexF);
68 
69  vout.general(m_vl, " Action_F_Ratio: %s\n", m_label.c_str());
70 
71  Field xi(NinF, NvolF, NexF);
72  rand->gauss_lex_global(xi);
73 
74  m_fopr_prec->set_config(m_U);
75  m_fopr->set_config(m_U);
76 
77  Field v1(NinF, NvolF, NexF), v2(NinF, NvolF, NexF);
78 
79  /*
80  m_fopr->set_mode("H");
81  v1 = m_fopr->mult(xi);
82 
83  int Nconv;
84  double diff;
85  m_fopr_prec->set_mode("DdagD");
86  m_solver_H_prec->solve(v2,v1,Nconv,diff);
87  vout.general(m_vl, " Nconv = %d diff = %.8e\n",Nconv,diff);
88 
89  m_fopr_prec->set_mode("H");
90  m_psf = m_fopr_prec->mult(v2);
91  */
92 
93  m_fopr->set_mode("H");
94  v2 = m_fopr->mult_dag(xi);
95 
96  m_fopr_prec->set_mode("H");
97  v1 = m_fopr_prec->mult_dag(v2);
98 
99  int Nconv;
100  double diff;
101  m_fopr_prec->set_mode("DdagD");
102  m_solver_H_prec->solve(m_psf, v1, Nconv, diff);
103  vout.general(m_vl, " Nconv = %d diff = %.8e\n", Nconv, diff);
104 
105  double xi2 = xi.norm();
106  double H_psf = xi2 * xi2;
107 
108  vout.general(m_vl, " H_Fratio = %18.8f\n", H_psf);
109  vout.general(m_vl, " H_F/dof = %18.8f\n", H_psf / size_psf);
110 
111  return H_psf;
112 }
113 
114 
115 //====================================================================
116 double Action_F_Ratio::calcH()
117 {
118  int Nvol = CommonParameters::Nvol();
119  int Ndim = CommonParameters::Ndim();
120 
121  int NinF = m_fopr_prec->field_nin();
122  int NvolF = m_fopr_prec->field_nvol();
123  int NexF = m_fopr_prec->field_nex();
124  int size_psf = NinF * NvolF * NexF * CommonParameters::NPE();
125 
126  Field v1(NinF, NvolF, NexF), v2(NinF, NvolF, NexF);
127 
128  vout.general(m_vl, " Action_F_Ratio: %s\n", m_label.c_str());
129 
130  m_fopr_prec->set_config(m_U);
131  m_fopr->set_config(m_U);
132 
133  m_fopr_prec->set_mode("H");
134  v1 = m_fopr_prec->mult(m_psf);
135 
136  int Nconv;
137  double diff;
138  m_fopr->set_mode("DdagD");
139  m_solver_H->solve(v2, v1, Nconv, diff);
140  vout.general(m_vl, " Nconv = %d diff = %.8e\n", Nconv, diff);
141 
142  double H_psf = v1 * v2;
143 
144  vout.general(m_vl, " H_Fratio = %18.8f\n", H_psf);
145  vout.general(m_vl, " H_F/dof = %18.8f\n", H_psf / size_psf);
146 
147  return H_psf;
148 }
149 
150 
151 //====================================================================
152 const Field Action_F_Ratio::force()
153 {
154  if (m_status_linkv == 0) {
155  int Nin = m_U->nin();
156  int Nvol = m_U->nvol();
157  int Nex = m_U->nex();
158  int Nc = CommonParameters::Nc();
159  int Ndim = CommonParameters::Ndim();
160 
161  int NinF = m_fopr_prec->field_nin();
162  int NvolF = m_fopr_prec->field_nvol();
163  int NexF = m_fopr_prec->field_nex();
164  Field eta(NinF, NvolF, NexF);
165 
166  vout.general(m_vl, " Action_F_Ratio: %s\n", m_label.c_str());
167 
168  m_fopr_prec->set_config(m_U);
169  m_fopr->set_config(m_U);
170  m_fopr_prec_force->set_config(m_U);
171  m_fopr_force->set_config(m_U);
172 
173  Field v1(NinF, NvolF, NexF), v2(NinF, NvolF, NexF);
174 
175  m_fopr_prec->set_mode("H");
176  v1 = m_fopr_prec->mult(m_psf);
177 
178  int Nconv;
179  double diff;
180  m_fopr->set_mode("DdagD");
181  m_solver_MD->solve(v2, v1, Nconv, diff);
182  vout.general(m_vl, " Solver: Nconv = %6d diff = %12.6e\n", Nconv, diff);
183 
184  Field force(Nin, Nvol, Nex);
185  Field force1(Nin, Nvol, Nex);
186 
187  force = m_fopr_force->force_core(v2);
188 
189  m_fopr_prec_force->set_mode("H");
190  force -= m_fopr_prec_force->force_core1(v2, m_psf);
191 
192  m_fopr_prec_force->set_mode("Hdag");
193  force -= m_fopr_prec_force->force_core1(m_psf, v2);
194 
195  m_force = force;
196  ++m_status_linkv;
197 
198  double Fave, Fmax, Fdev;
199  m_force.stat(Fave, Fmax, Fdev);
200  vout.general(m_vl, " Fratio_ave = %12.6f Fratio_max = %12.6f Fratio_dev = %12.6f\n",
201  Fave, Fmax, Fdev);
202 
203  return m_force;
204  } else {
205  vout.general(m_vl, " Action_F_Ratio returns previous force.\n");
206  return m_force;
207  }
208 }
209 
210 
211 //====================================================================
212 //============================================================END=====