Bridge++  Ver. 1.2.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
projection_Maximum_SU_N.cpp
Go to the documentation of this file.
1 
15 
16 #ifdef USE_PARAMETERS_FACTORY
17 #include "parameters_factory.h"
18 #endif
19 
20 #ifdef USE_FACTORY
21 namespace {
22  Projection *create_object()
23  {
24  return new Projection_Maximum_SU_N();
25  }
26 
27 
28  bool init = Projection::Factory::Register("Maximum_SU_N", create_object);
29 }
30 #endif
31 
32 //- parameter entries
33 namespace {
34  void append_entry(Parameters& param)
35  {
36  param.Register_int("maximum_number_of_iteration", 0);
37  param.Register_double("convergence_criterion", 0.0);
38 
39  param.Register_string("verbose_level", "NULL");
40  }
41 
42 
43 #ifdef USE_PARAMETERS_FACTORY
44  bool init_param = ParametersFactory::Register("Projection.Maximum_SU_N", append_entry);
45 #endif
46 }
47 //- end
48 
49 //- parameters class
51 //- end
52 
53 const std::string Projection_Maximum_SU_N::class_name = "Projection_Maximum_SU_N";
54 
55 //====================================================================
57 {
58  const string str_vlevel = params.get_string("verbose_level");
59 
60  m_vl = vout.set_verbose_level(str_vlevel);
61 
62  //- fetch and check input parameters
63  int Niter;
64  double Enorm;
65 
66  int err = 0;
67  err += params.fetch_int("maximum_number_of_iteration", Niter);
68  err += params.fetch_double("convergence_criterion", Enorm);
69 
70  if (err) {
71  vout.crucial(m_vl, "%s: fetch error, input parameter not found.\n", class_name.c_str());
72  abort();
73  }
74 
75 
76  set_parameters(Niter, Enorm);
77 }
78 
79 
80 //====================================================================
81 void Projection_Maximum_SU_N::set_parameters(const int Niter, const double Enorm)
82 {
83  //- print input parameters
84  vout.general(m_vl, "Parameters of %s:\n", class_name.c_str());
85  vout.general(m_vl, " Niter = %d\n", Niter);
86  vout.general(m_vl, " Enorm = %12.4e\n", Enorm);
87 
88  //- range check
89  int err = 0;
90  err += ParameterCheck::non_negative(Niter);
91  err += ParameterCheck::non_zero(Enorm);
92 
93  if (err) {
94  vout.crucial(m_vl, "%s: parameter range check failed.\n", class_name.c_str());
95  abort();
96  }
97 
98  //- store values
99  m_Niter = Niter;
100  m_Enorm = Enorm;
101 }
102 
103 
104 //====================================================================
106  double alpha,
107  const Field_G& Cst, const Field_G& Uorg)
108 {
109  int Nex = Uorg.nex();
110  int Nvol = Uorg.nvol();
111  int Nc = CommonParameters::Nc();
112 
113  assert(Cst.nex() == Nex);
114  assert(Cst.nvol() == Nvol);
115  assert(U.nex() == Nex);
116  assert(U.nvol() == Nvol);
117 
118  Field_G u_tmp(Nvol, Nex);
119  for (int ex = 0; ex < Nex; ++ex) {
120  u_tmp.setpart_ex(ex, Cst, ex);
121  u_tmp.addpart_ex(ex, Uorg, ex, 1.0 - alpha);
122  }
123 
124  maxTr(U, u_tmp);
125 }
126 
127 
128 //====================================================================
130  double alpha, const Field_G& Sigmap,
131  const Field_G& Cst, const Field_G& Uorg)
132 {
133  vout.crucial(m_vl, "%s: force_recursive() is not available.\n", class_name.c_str());
134  abort();
135 }
136 
137 
138 //====================================================================
140 {
141  int Nc = CommonParameters::Nc();
142  int Nvol = Cst.nvol();
143  int Nex = Cst.nex();
144 
145  assert(Nvol == G0.nvol());
146  assert(Nex == G0.nex());
147 
148  int Nmt = 1; // number of subgroup maximization loop:
149  // seems not important because of outer iter-loop.
150 
151  Mat_SU_N unity(Nc);
152  unity.unit();
153 
154  Field_G Udelta(Nvol, Nex), A(Nvol, Nex);
155  A = Cst;
156 
157  vout.detailed(m_vl, "Maximum projection start.\n");
158 
159  for (int ex = 0; ex < Nex; ++ex) {
160  for (int site = 0; site < Nvol; ++site) {
161  G0.set_mat(site, ex, unity);
162  }
163  }
164 
165  for (int iter = 0; iter < m_Niter; ++iter) {
166  for (int ex = 0; ex < Nex; ++ex) {
167  for (int site = 0; site < Nvol; ++site) {
168  Udelta.set_mat(site, ex, unity);
169  }
170  }
171 
172  for (int imt = 0; imt < Nmt; ++imt) {
173  for (int i1 = 0; i1 < Nc; ++i1) {
174  int i2 = (i1 + 1) % Nc;
175  maxTr_SU2(i1, i2, G0, A, Udelta);
176  }
177  }
178  // for exact check with Fortran version (passed).
179 
180  /*
181  for(int imt = 0; imt < Nmt; ++imt){
182  for(int i = Nc; i > 0; --i){
183  int i1 = i % Nc;
184  int i2 = (i1 + 1) % Nc;
185  maxTr_SU2(i1, i2, G0, A, Udelta);
186  }
187  }
188  */
189 
190  //- convergence test
191  double retr1 = 0.0;
192  for (int ex = 0; ex < Nex; ++ex) {
193  for (int site = 0; site < Nvol; ++site) {
194  for (int cc = 0; cc < Nc; ++cc) {
195  retr1 += Udelta.cmp_r(cc * (1 + Nc), site, ex);
196  }
197  }
198  }
199 
200  double retr = Communicator::reduce_sum(retr1);
201  int Npe = Communicator::size();
202  double deltaV = 1.0 - retr / (Nc * Nvol * Nex * Npe);
203  vout.detailed(m_vl, " iter = %d deltaV = %12.4e\n", iter, deltaV);
204 
205  if (deltaV < m_Enorm) {
206  Mat_SU_N ut(Nc);
207  for (int ex = 0; ex < Nex; ++ex) {
208  for (int site = 0; site < Nvol; ++site) {
209  G0.mat_dag(ut, site, ex);
210  G0.set_mat(site, ex, ut);
211  }
212  }
213 
214  vout.detailed(m_vl, "Maximum projection converged.\n");
215 
216  return;
217  }
218  }
219 
220 
221  vout.crucial(m_vl, "Maximum projection not converged.\n");
222  abort();
223 }
224 
225 
226 //====================================================================
227 void Projection_Maximum_SU_N::maxTr_SU2(int i1, int i2, Field_G& Gmax,
228  Field_G& A, Field_G& Udelta)
229 {
230  int Nc = CommonParameters::Nc();
231  int Nvol = A.nvol();
232  int Nex = A.nex();
233 
234  assert(i1 < Nc);
235  assert(i2 < Nc);
236 
237  int j1 = mindex(i1, i1, Nc);
238  int j2 = mindex(i2, i2, Nc);
239  int k1 = mindex(i1, i2, Nc);
240  int k2 = mindex(i2, i1, Nc);
241 
242  Mat_SU_N at(Nc), vt(Nc);
243  Field_G v(Nvol, Nex), w(Nvol, Nex);
244 
245  //----------[ | # # 0 | <i1 ]--------------------------
246  //----------[ V = | # # 0 | <i2 ]--------------------------
247  //----------[ | 0 0 1 | ]--------------------------
248 
249  for (int ex = 0; ex < Nex; ++ex) {
250  for (int site = 0; site < Nvol; ++site) {
251  at = A.mat(site, ex);
252 
253  double xlamd =
254  at.r(j1) * at.r(j1) + at.i(j1) * at.i(j1) + 2.0 * at.r(j1) * at.r(j2)
255  + at.r(k1) * at.r(k1) + at.i(k1) * at.i(k1) - 2.0 * at.i(j1) * at.i(j2)
256  + at.r(k2) * at.r(k2) + at.i(k2) * at.i(k2) - 2.0 * at.r(k1) * at.r(k2)
257  + at.r(j2) * at.r(j2) + at.i(j2) * at.i(j2) + 2.0 * at.i(k1) * at.i(k2);
258  xlamd = 1.0 / sqrt(xlamd);
259 
260  vt.unit();
261  vt.set(j1, (at.r(j1) + at.r(j2)) * xlamd, (-at.i(j1) + at.i(j2)) * xlamd);
262  vt.set(k1, (at.r(k2) - at.r(k1)) * xlamd, (-at.i(k2) - at.i(k1)) * xlamd);
263  vt.set(k2, (at.r(k1) - at.r(k2)) * xlamd, (-at.i(k1) - at.i(k2)) * xlamd);
264  vt.set(j2, (at.r(j1) + at.r(j2)) * xlamd, (at.i(j1) - at.i(j2)) * xlamd);
265 
266  v.set_mat(site, ex, vt);
267  }
268  }
269 
270  for (int ex = 0; ex < Nex; ++ex) {
271  mult_Field_Gnn(w, ex, A, ex, v, ex);
272  }
273  A = w;
274 
275  for (int ex = 0; ex < Nex; ++ex) {
276  mult_Field_Gnn(w, ex, Gmax, ex, v, ex);
277  }
278  Gmax = w;
279 
280  for (int ex = 0; ex < Nex; ++ex) {
281  mult_Field_Gnn(w, ex, Udelta, ex, v, ex);
282  }
283  Udelta = w;
284 }
285 
286 
287 //====================================================================
288 //============================================================END=====
BridgeIO vout
Definition: bridgeIO.cpp:207
void detailed(const char *format,...)
Definition: bridgeIO.cpp:50
void Register_string(const string &, const string &)
Definition: parameters.cpp:352
void general(const char *format,...)
Definition: bridgeIO.cpp:38
void Register_int(const string &, const int)
Definition: parameters.cpp:331
void set_parameters(const Parameters &params)
void maxTr(Field_G &U, Field_G &V)
double m_Enorm
convergence criterion of maximization
int nvol() const
Definition: field.h:101
Class for parameters.
Definition: parameters.h:40
SU(N) gauge field.
Definition: field_G.h:36
void mult_Field_Gnn(Field_G &w, const int ex, const Field_G &u1, const int ex1, const Field_G &u2, const int ex2)
int nex() const
Definition: field.h:102
void project(Field_G &U, double alpha, const Field_G &C, const Field_G &Uorg)
projection U = P[alpha, C, Uorg]
int m_Niter
maximum iteration of maximization steps
void crucial(const char *format,...)
Definition: bridgeIO.cpp:26
base class for projection operator into gauge group.
Definition: projection.h:33
static const std::string class_name
static bool Register(const std::string &realm, const creator_callback &cb)
static int size()
size of small world.
void maxTr_SU2(int, int, Field_G &, Field_G &, Field_G &)
int mindex(int i, int j, int Nc)
int non_zero(const double v)
Definition: checker.cpp:31
void force_recursive(Field_G &Xi, Field_G &iTheta, double alpha, const Field_G &Sigmap, const Field_G &C, const Field_G &U)
force calculation: invalid in this class.
Mat_SU_N mat_dag(const int site, const int mn=0) const
Definition: field_G.h:123
int non_negative(const int v)
Definition: checker.cpp:21
Mat_SU_N & unit()
Definition: mat_SU_N.h:373
static int reduce_sum(int count, double *recv_buf, double *send_buf, int pattern=0)
make a global sum of an array of double over the communicator. pattern specifies the dimensions to be...
void Register_double(const string &, const double)
Definition: parameters.cpp:324
double r(int c) const
Definition: mat_SU_N.h:114
void set(int c, double re, const double &im)
Definition: mat_SU_N.h:133
int fetch_double(const string &key, double &val) const
Definition: parameters.cpp:124
string get_string(const string &key) const
Definition: parameters.cpp:85
void set_mat(const int site, const int mn, const Mat_SU_N &U)
Definition: field_G.h:156
Mat_SU_N mat(const int site, const int mn=0) const
Definition: field_G.h:110
int fetch_int(const string &key, int &val) const
Definition: parameters.cpp:141
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:191
Bridge::VerboseLevel m_vl
Definition: projection.h:36