Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
test_Spectrum_CRSMatrix_Overlap5d.cpp
Go to the documentation of this file.
1 
14 #include "parameterManager_YAML.h"
15 
16 #include "bridgeIO.h"
17 using Bridge::vout;
18 
19 // #include "index_lex.h"
20 #include "index_eo.h"
21 
22 #include "fopr_CRS.h"
23 #include "fopr_Wilson_eo.h"
24 #include "fopr_Overlap_5d.h"
25 #include "fopr_Wilson.h"
26 #include "fopr_Overlap.h"
27 
28 #include "gaugeConfig.h"
29 
30 #include "source_4spinor_Local.h"
31 
32 #include "solver_CG.h"
33 
34 #ifdef USE_TEST
35 #include "test.h"
36 #endif
37 
38 #ifdef USE_TESTMANAGER_AUTOREGISTER
39 #include "testManager.h"
40 #endif
41 
42 //====================================================================
44 
59 namespace Test_Spectrum_CRSMatrix {
60  //- test-private parameters
61  namespace {
62  const std::string filename_input = "test_Spectrum_CRSMatrix_Overlap5d.yaml";
63  const std::string filename_output = "stdout";
64 
65  class Parameters_Test_Spectrum_CRSMatrix : public Parameters {
66  public:
67  Parameters_Test_Spectrum_CRSMatrix()
68  {
69  Register_string("gauge_config_type_input", "NULL");
70  Register_string("config_filename_input", "NULL");
71 
72  Register_string("matrix_output", "NULL");
73  Register_string("source_output", "NULL");
74  Register_string("solution_output", "NULL");
75 
76  Register_string("verbose_level", "NULL");
77 
78  Register_double("expected_result", 0.0);
79  }
80  };
81  }
82 
83  //- prototype declaration
84  int overlap_5d(void);
85 
86  int CRSsolver(
87  const string& solution,
88  const string& matrix,
89  const string& source,
90  double& result /* return value */
91  );
92 
93 
94 #ifdef USE_TESTMANAGER_AUTOREGISTER
95 #ifdef USE_MPI
96  // this test runs only in single-node environment.
97 #else
98  //- NB. CRS test is skipped, because it is time-consuming.
99  // namespace {
100  // static const bool is_registered = TestManager::RegisterTest(
101  // "Spectrum.CRSMatrix.Overlap5d",
102  // overlap_5d
103  // );
104  // }
105 #endif
106 #endif
107 
108  //====================================================================
109  int overlap_5d(void)
110  {
111  // #### parameter setup ####
112  int Nc = CommonParameters::Nc();
113  int Nd = CommonParameters::Nd();
114  int Nvol = CommonParameters::Nvol();
115  int Ndim = CommonParameters::Ndim();
116 
117  Parameters_Test_Spectrum_CRSMatrix params_test;
118 
119  Parameters params_all;
120 
121  params_all.Register_Parameters("Test_Spectrum_CRSMatrix", &params_test);
122 
123  ParameterManager_YAML params_manager;
124  params_manager.read_params(filename_input, &params_all);
125 
126  const string str_gconf_read = params_test.get_string("gauge_config_type_input");
127  const string readfile = params_test.get_string("config_filename_input");
128  const string matrix_file = params_test.get_string("matrix_output");
129  const string source_file = params_test.get_string("source_output");
130  const string solution_file = params_test.get_string("solution_output");
131  const string str_vlevel = params_test.get_string("verbose_level");
132 #ifdef USE_TEST
133  const double expected_result = params_test.get_double("expected_result");
134 #endif
135 
136  Bridge::VerboseLevel vl = vout.set_verbose_level(str_vlevel);
137 
138  //- print input parameters
139  vout.general(vl, " gconf_read = %s\n", str_gconf_read.c_str());
140  vout.general(vl, " readfile = %s\n", readfile.c_str());
141  vout.general(vl, " matrix_output = %s\n", matrix_file.c_str());
142  vout.general(vl, " source_output = %s\n", source_file.c_str());
143  vout.general(vl, " solution_output = %s\n", solution_file.c_str());
144  vout.general(vl, " vlevel = %s\n", str_vlevel.c_str());
145 
146  //- input parameter check
147  int err = 0;
148  err += ParameterCheck::non_NULL(str_gconf_read);
149  err += ParameterCheck::non_NULL(readfile);
150  err += ParameterCheck::non_NULL(matrix_file);
151  err += ParameterCheck::non_NULL(source_file);
152  err += ParameterCheck::non_NULL(solution_file);
153 
154  if (err) {
155  vout.crucial(vl, "Test_Spectrum_CRSMatrix: Input parameters have not been set.\n");
156  abort();
157  }
158 
159 
160  // #### Set up a gauge configuration ####
161  Field_G *U = new Field_G(Nvol, Ndim);
162  GaugeConfig *gconf_read = new GaugeConfig(str_gconf_read);
163  gconf_read->read_file((Field *)U, readfile);
164  // gconf->set_cold((Field*)U);
165 
166  Field_G *Ueo = new Field_G(Nvol, Ndim);
167  Index_eo index_eo;
168  index_eo.convertField(*Ueo, *U);
169 
170 
171  // #### object setup #####
172  //- overlap parameters
173  double mq = 0.2;
174  double M0 = 1.6;
175  int Np = 4;
176  double x_min = 0.10;
177  double x_max = 8.0;
178  double Niter_ms = 1000;
179  double Stop_cond_ms = 1.0e-28;
180  valarray<int> boundary(Ndim);
181  boundary[0] = 1;
182  boundary[1] = 1;
183  boundary[2] = 1;
184  boundary[3] = -1;
185 
186 
187  Fopr_Wilson_eo *fopr_w_eo = new Fopr_Wilson_eo("Dirac");
188  double kappa = 0.5 / (4.0 - M0);
189  fopr_w_eo->set_parameters(kappa, boundary);
190  fopr_w_eo->set_config(Ueo);
191 
192  Fopr_Overlap_5d *fopr_ov5d = new Fopr_Overlap_5d(fopr_w_eo);
193  fopr_ov5d->set_parameters(mq, M0, Np, x_min, x_max,
194  Niter_ms, Stop_cond_ms, boundary);
195  fopr_ov5d->set_mode("D_eo");
196 
197  Fopr_CRS *fopr_crs = new Fopr_CRS(fopr_ov5d);
198  fopr_crs->write_matrix(matrix_file);
199 
200 
201  //- NB. fopr_ov is used only for check
202  Fopr_Wilson *fopr_w = new Fopr_Wilson;
203  fopr_w->set_parameters(kappa, boundary);
204  fopr_w->set_config(U);
205 
206  Fopr_Overlap *fopr_ov = new Fopr_Overlap(fopr_w);
207  fopr_ov->set_parameters(mq, M0, Np, x_min, x_max,
208  Niter_ms, Stop_cond_ms, boundary);
209 
210 
211  int Niter_5d = 1000;
212  double Stop_cond_5d = 1.0e-28;
213  Solver *solver = new Solver_CG(fopr_ov5d);
214  solver->set_parameters(Niter_5d, Stop_cond_5d);
215 
216  valarray<int> source_position(Ndim);
217  source_position[0] = 0;
218  source_position[1] = 0;
219  source_position[2] = 0;
220  source_position[3] = 0;
221  Source_4spinor_Local *source = new Source_4spinor_Local();
222  source->set_parameters(source_position);
223 
224 
225  // #### Execution main part ####
226  valarray<Field_F> sq(Nc * Nd);
227  for (int cd = 0; cd < Nc * Nd; ++cd) {
228  sq[cd] = 0.0;
229  }
230 
231  Field_F b;
232  Field xq((Field)b);
233  Field v1(xq);
234 
235  int Nin = fopr_ov5d->field_nin();
236  int Nvol2 = fopr_ov5d->field_nvol();
237  int Npl = fopr_ov5d->field_nex();
238 
239  Field v(Nin, Nvol2, Npl);
240  Field t(Nin, Nvol2, Npl);
241  Field p(Nin, Nvol2, Npl);
242  Field s(Nin, Nvol2, Npl);
243  Field x(Nin, Nvol2, Npl);
244 
245  Field be(Nin, Nvol2, 1);
246  Field bo(Nin, Nvol2, 1);
247  Field xqe(Nin, Nvol2, 1);
248  Field xqo(Nin, Nvol2, 1);
249  Field xt(Nin, Nvol2, 1);
250 
251  int Nconv;
252  double diff;
253 
254  double vv = 0.0L;
255 
256  // { int ispin = 0;
257  // { int icolor = 0;
258  for (int ispin = 0; ispin < Nd; ++ispin) {
259  for (int icolor = 0; icolor < Nc; ++icolor) {
260  source->set(b, icolor, ispin);
261 
262  double snorm = 1.0 / b.norm2();
263 
264  index_eo.convertField(be, b, 0);
265  index_eo.convertField(bo, b, 1);
266 
267  v = 0.0;
268  v.setpart_ex(2 * Np, be, 0);
269  fopr_ov5d->LUprecond(t, v, 0);
270 
271  v = 0.0;
272  v.setpart_ex(2 * Np, bo, 0);
273  fopr_ov5d->LUprecond(p, v, 1);
274  fopr_ov5d->Mopr_5d_eo(v, p, 0);
275  fopr_ov5d->LUprecond(p, v, 0);
276 
277  t -= p;
278  t.write_text(source_file);
279  s = fopr_ov5d->DD_5d_eo(t, -1);
280 
281 
282  fopr_ov5d->set_mode("DdagD_eo");
283  solver->solve(x, s, Nconv, diff);
284 
285  x.write_text(solution_file);
286 
287  //- Solution
288  xqe.setpart_ex(0, x, 2 * Np);
289 
290  v = 0.0;
291  v.setpart_ex(2 * Np, bo, 0);
292  fopr_ov5d->LUprecond(p, v, 1);
293  xqo.setpart_ex(0, p, 2 * Np);
294 
295  fopr_ov5d->Mopr_5d_eo(p, x, 1);
296  fopr_ov5d->LUprecond(t, p, 1);
297  xt.setpart_ex(0, t, 2 * Np);
298  xqo -= xt;
299 
300  index_eo.reverseField(xq, xqe, 0);
301  index_eo.reverseField(xq, xqo, 1);
302 
303  fopr_ov->set_mode("H");
304  v1 = fopr_ov->mult(xq);
305  v1 -= (Field)b;
306  vv = v1 * v1;
307  vv = vv * snorm;
308  vout.general(vl, " Diff(ov) = %16.8e\n", vv);
309 
310  sq[icolor + Nc * ispin] = (Field_F)xq;
311  }
312  }
313 
314  double result;
315 
316  CRSsolver(solution_file, matrix_file, source_file, result);
317 
318 
319  // #### tydy up ####
320  delete U;
321  delete gconf_read;
322  delete Ueo;
323 
324  delete fopr_crs;
325  delete fopr_ov;
326  delete fopr_w;
327 
328  delete fopr_ov5d;
329  delete fopr_w_eo;
330 
331  delete solver;
332  delete source;
333 
334 
335 #ifdef USE_TEST
336  return Test::verify(expected_result, result);
337 
338 #else
339  return EXIT_SUCCESS;
340 #endif
341  }
342 } // namespace Test_Spectrum_CRSMatrix