CHROMA
symm_prec_tests.cc
Go to the documentation of this file.
1 #include "chromabase.h"
2 
3 #include "handle.h"
4 
5 #include "seoprec_linop.h"
6 #include "eoprec_linop.h"
7 #include "eoprec_wilstype_fermact_w.h"
8 #include "seoprec_wilstype_fermact_w.h"
9 #include "actions/ferm/linop/eoprec_clover_linop_w.h"
10 #include "actions/ferm/linop/seoprec_clover_linop_w.h"
11 #include "actions/ferm/fermacts/fermact_factory_w.h"
12 #include "util/gauge/reunit.h"
13 #include "gtest/gtest.h"
14 
15 #include "io/xml_group_reader.h"
16 #include "./symm_prec_xml.h"
17 
18 #ifdef BUILD_QUDA
19 #include "quda.h"
20 #include "update/molecdyn/predictor/quda_predictor.h"
21 #endif
22 
23 using namespace Chroma;
24 using namespace QDP;
25 using namespace SymmPrecTesting;
26 
27 template<typename TestType>
28 class SymmFixtureT : public TestType {
29 public:
30  using T = LatticeFermion;
31  using Q = multi1d<LatticeColorMatrix>;
32  using P = multi1d<LatticeColorMatrix>;
33 
34  using S_asymm_T = EvenOddPrecWilsonTypeFermAct<T,P,Q>;
35  using S_symm_T = SymEvenOddPrecWilsonTypeFermAct<T,P,Q>;
36  using LinOpSymm_T = SymEvenOddPrecCloverLinOp;
37  using LinOpAsymm_T = EvenOddPrecCloverLinOp;
38 
39 
40  void SetUp() {
41  u.resize(Nd);
42  for(int mu=0; mu < Nd; ++mu) {
43  gaussian(u[mu]);
44  reunit(u[mu]);
45  }
46 
47  std::istringstream input_asymm(fermact_xml_asymm);
48  std::istringstream input_symm(fermact_xml_symm);
49 
50  XMLReader xml_in_asymm(input_asymm);
51  XMLReader xml_in_symm(input_symm);
52 
53  std::istringstream input_asymm_periodic(fermact_xml_asymm_periodic);
54  std::istringstream input_symm_periodic(fermact_xml_symm_periodic);
55  std::istringstream twisted_fermact(fermact_xml_symm_twisted);
56 
57  XMLReader xml_in_asymm_periodic(input_asymm_periodic);
58  XMLReader xml_in_symm_periodic(input_symm_periodic);
59  XMLReader xml_in_twisted(twisted_fermact);
60 
61  S_asymm = dynamic_cast<S_asymm_T*>(TheFermionActionFactory::Instance().createObject("CLOVER",
62  xml_in_asymm,
63  "FermionAction"));
64 
65  S_symm = dynamic_cast<S_symm_T*>(TheFermionActionFactory::Instance().createObject("SEOPREC_CLOVER",
66  xml_in_symm,
67  "FermionAction"));
68 
69  S_asymm_periodic = dynamic_cast<S_asymm_T*>(TheFermionActionFactory::Instance().createObject("CLOVER",
70  xml_in_asymm_periodic,
71  "FermionAction"));
72 
73  S_symm_periodic = dynamic_cast<S_symm_T*>(TheFermionActionFactory::Instance().createObject("SEOPREC_CLOVER",
74  xml_in_symm_periodic,
75  "FermionAction"));
76 
77 
78  S_symm_twisted = dynamic_cast<S_symm_T*>(TheFermionActionFactory::Instance().createObject("SEOPREC_CLOVER",
79  xml_in_twisted,
80  "FermionAction"));
81  state = S_asymm->createState(u);
82 
83  M_asymm =dynamic_cast<LinOpAsymm_T *>(S_asymm->linOp(state));
84  M_symm =dynamic_cast<LinOpSymm_T *>(S_symm->linOp(state));
85  M_tw = dynamic_cast<LinOpSymm_T *>(S_symm_twisted->linOp(state));
86  }
87 
88 
89  void TearDown() {}
90 
91  Q u;
92  Handle<S_symm_T> S_symm;
93  Handle<S_asymm_T> S_asymm;
94  Handle<S_symm_T> S_symm_periodic;
95  Handle<S_asymm_T> S_asymm_periodic;
96  Handle<S_symm_T> S_symm_twisted;
97  Handle<FermState<T,P,Q> > state;
98 
99  Handle<LinOpAsymm_T> M_asymm;
100  Handle<LinOpSymm_T> M_symm;
101  Handle<LinOpSymm_T> M_tw;
102 };
103 
104 class SymmFixture : public SymmFixtureT<::testing::Test> {};
105 class QPropTest : public SymmFixtureT<::testing::TestWithParam<std::string>>{};
106 
107 TEST_F(SymmFixture, CheckOp)
108 {
109  // Verify A_oo M_symm = M_asymm
110  LatticeFermion x=zero;
111  gaussian(x, rb[1]);
112 
113  LatticeFermion t1 = zero;
114 
115 
116  LatticeFermion AM_symm_x = zero;
117  LatticeFermion M_asymm_x = zero;
118 
119  {
120  QDPIO::cout << "Op: " << std::endl;
121  // AM_symm_x = A_oo t_1 = A_oo M_symm x
122  (*M_symm)(t1,x,PLUS);
123  (*M_symm).unprecOddOddLinOp(AM_symm_x,t1,PLUS);
124 
125 
126  (*M_asymm)(M_asymm_x, x, PLUS);
127 
128  M_asymm_x -= AM_symm_x;
129  Double norm_cb0 = sqrt(norm2(M_asymm_x,rb[0]));
130  Double norm_cb1 = sqrt(norm2(M_asymm_x,rb[1]));
131 
132  QDPIO::cout << "CB=0: || A_oo M_symm_x - M_asymm_x ||= " << norm_cb0 <<std::endl;
133  QDPIO::cout << "CB=1: || A_oo M_symm_x - M_asymm_x ||= " << norm_cb1
134  << " || A_oo M_symm_x - M_asymm_x ||/ ||x||=" << norm_cb1/sqrt(norm2(x,rb[1]))
135  << std::endl;
136 
137  ASSERT_LT( toDouble(norm_cb0), 1.0e-14);
138  ASSERT_LT( toDouble(norm_cb1), 1.0e-13);
139  }
140  // Now check herm conj
141  {
142  QDPIO::cout << "Daggered Op: " << std::endl;
143  (*M_symm).unprecOddOddLinOp(t1,x,MINUS);
144  (*M_symm)(AM_symm_x,t1,MINUS);
145 
146  (*M_asymm)(M_asymm_x, x, MINUS);
147  M_asymm_x -= AM_symm_x;
148  Double norm_cb0 = sqrt(norm2(M_asymm_x,rb[0]));
149  Double norm_cb1 = sqrt(norm2(M_asymm_x,rb[1]));
150 
151  QDPIO::cout << "CB=0: || A_oo M_symm_x - M_asymm_x ||= " << norm_cb0 <<std::endl;
152  QDPIO::cout << "CB=1: || A_oo M_symm_x - M_asymm_x ||= " << norm_cb1
153  << " || A_oo M_symm_x - M_asymm_x ||/ ||x||=" << norm_cb1/sqrt(norm2(x,rb[1]))
154  << std::endl;
155 
156  ASSERT_LT( toDouble(norm_cb0), 1.0e-14);
157  ASSERT_LT( toDouble(norm_cb1), 1.0e-13);
158  }
159 
160 }
161 
162 // Check both symm and asymm linops have same unprec op
163 TEST_F(SymmFixture, CheckUnprecOp)
164 {
165  T x;
166  T unprec_symm_x;
167  T unprec_asymm_x;
168 
169  gaussian(x);
170  {
171 
172  QDPIO::cout << "Regular Op:" << std::endl;
173 
174  (*M_symm).unprecLinOp( unprec_symm_x, x, PLUS);
175  (*M_asymm).unprecLinOp( unprec_asymm_x, x, PLUS);
176 
177  unprec_asymm_x -= unprec_symm_x;
178  Double norm_cb0 = sqrt(norm2(unprec_asymm_x,rb[0]));
179  Double norm_cb1 = sqrt(norm2(unprec_asymm_x,rb[1]));
180 
181  QDPIO::cout << "CB=0: || Munprec_symm_x - Munprec_asymm_x ||= " << norm_cb0 <<std::endl;
182  QDPIO::cout << "CB=1: || Munprec_symm_x - Munprec_asymm_x ||= " << norm_cb1 << std::endl;
183 
184  ASSERT_LT( toDouble(norm_cb0), 1.0e-13);
185  ASSERT_LT( toDouble(norm_cb1), 1.0e-13);
186 
187  }
188 
189  {
190  QDPIO::cout << "Daggered Op:" << std::endl;
191 
192  (*M_symm).unprecLinOp( unprec_symm_x, x, MINUS);
193  (*M_asymm).unprecLinOp( unprec_asymm_x, x, MINUS);
194  unprec_asymm_x -= unprec_symm_x;
195 
196  Double norm_cb0 = sqrt(norm2(unprec_asymm_x,rb[0]));
197  Double norm_cb1 = sqrt(norm2(unprec_asymm_x,rb[1]));
198 
199  QDPIO::cout << "CB=0: || Munprec_symm_x - Munprec_asymm_x ||= " << norm_cb0 <<std::endl;
200  QDPIO::cout << "CB=1: || Munprec_symm_x - Munprec_asymm_x ||= " << norm_cb1 << std::endl;
201 
202  ASSERT_LT( toDouble(norm_cb0), 1.0e-13);
203  ASSERT_LT( toDouble(norm_cb1), 1.0e-13);
204 
205  }
206 
207 }
208 
209 // Check QProp Functionality.
210 TEST_P(QPropTest, CheckQprop)
211 {
212  LatticeFermion rhs=zero;
213  gaussian(rhs);
214 
215  std::istringstream inv_param_xml_stream(GetParam());
216  XMLReader xml_in(inv_param_xml_stream);
217 
218  GroupXML_t inv_param = readXMLGroup(xml_in, "//InvertParam", "invType");
219  Handle<SystemSolver<T>> qprop_solver = S_symm->qprop(state,inv_param);
220  LatticeFermion x = zero;
221 
222  (*qprop_solver)(x,rhs);
223 
224  // Check residuum
225  LatticeFermion Ax=zero;
226  (*M_symm).unprecLinOp(Ax,x,PLUS);
227  Ax -= rhs;
228 
229  Double resid_cb0 = sqrt(norm2(Ax,rb[0]));
230  Double resid_cb1 = sqrt(norm2(Ax,rb[1]));
231  QDPIO::cout << "Qprop: rsd cb0 = " << resid_cb0 << std::endl;
232  QDPIO::cout << "Qprop: rsd cb1 = " << resid_cb1 << std::endl;
233 
234  Double resid = sqrt(norm2(Ax));
235  Double resid_rel = resid/sqrt(norm2(rhs));
236  QDPIO::cout << "QProp Check Back: || r || = " << resid << " || r ||/||b|| = "
237  << resid_rel << std::endl;
238 
239  ASSERT_LT(toDouble(resid_rel), 1.0e-8);
240 
241 }
242 
243 INSTANTIATE_TEST_CASE_P(PropSyssolver,
244  QPropTest,
245  ::testing::Values(inv_param_syssolver_bicgstab_xml));
246 
247 #ifdef BUILD_QUDA
248 INSTANTIATE_TEST_CASE_P(PropQUDASolver,
249  QPropTest,
250  ::testing::Values(inv_param_quda_bicgstab_xml,
251  inv_param_quda_multigrid_xml));
252 
253 #endif
254 
255 class MdagMInvTestSymm : public SymmFixtureT<::testing::TestWithParam<std::string>>{};
256 class MdagMInvTestAsymm : public SymmFixtureT<::testing::TestWithParam<std::string>>{};
257 
258 TEST_P(MdagMInvTestSymm, CheckMdagMInvSymm)
259 {
260  std::istringstream inv_param_xml_stream(GetParam());
261  XMLReader xml_in(inv_param_xml_stream);
262 
263  GroupXML_t inv_param = readXMLGroup(xml_in, "//InvertParam", "invType");
264  Handle<MdagMSystemSolver<T>> MdagM_solver = S_symm->invMdagM(state,inv_param);
265 
266  T b; gaussian(b,rb[1]);
267  T x; x[rb[1]] = zero;
268 
269  (*MdagM_solver)(x,b);
270 
271  T tmp, r;
272 
273  (*M_symm)(tmp,x,PLUS);
274  (*M_symm)(r,tmp,MINUS);
275 
276  r[rb[1]] -= b;
277  Double resid = sqrt(norm2(r,rb[1]));
278  Double resid_rel = resid/sqrt(norm2(b,rb[1]));
279  QDPIO::cout << "MdagM check: || r || = " << resid << " || r || / || b ||=" << resid_rel << std::endl;
280  ASSERT_LT( toDouble(resid_rel),1.0e-8);
281 
282 
283 
284 }
285 
286 #ifdef BUILD_QUDA
287 TEST_P(MdagMInvTestSymm, CheckMdagMInvSymmQUDAPredict)
288 {
289  std::istringstream inv_param_xml_stream(GetParam());
290  XMLReader xml_in(inv_param_xml_stream);
291 
292  GroupXML_t inv_param = readXMLGroup(xml_in, "//InvertParam", "invType");
293  Handle<MdagMSystemSolver<T>> MdagM_solver = S_symm->invMdagM(state,inv_param);
294 
295  T b; gaussian(b,rb[1]);
296  T x; x[rb[1]] = zero;
297  QUDA4DChronoPredictor chrono(5,DEFAULT);
298  (*MdagM_solver)(x,b,chrono);
299 
300  T tmp, r;
301 
302  (*M_symm)(tmp,x,PLUS);
303  (*M_symm)(r,tmp,MINUS);
304 
305  r[rb[1]] -= b;
306  Double resid = sqrt(norm2(r,rb[1]));
307  Double resid_rel = resid/sqrt(norm2(b,rb[1]));
308  QDPIO::cout << "MdagM check: || r || = " << resid << " || r || / || b ||=" << resid_rel << std::endl;
309  ASSERT_LT( toDouble(resid_rel),1.0e-8);
310 
311 
312 
313 }
314 #endif
315 
316 TEST_P(MdagMInvTestAsymm, CheckMdagMInvAsymm)
317 {
318  std::istringstream inv_param_xml_stream(GetParam());
319  XMLReader xml_in(inv_param_xml_stream);
320 
321  GroupXML_t inv_param = readXMLGroup(xml_in, "//InvertParam", "invType");
322  Handle<MdagMSystemSolver<T>> MdagM_solver = S_asymm->invMdagM(state,inv_param);
323 
324  T b; gaussian(b,rb[1]);
325  T x; x[rb[1]] = zero;
326 
327  (*MdagM_solver)(x,b);
328 
329  T tmp, r;
330 
331  (*M_asymm)(tmp,x,PLUS);
332  (*M_asymm)(r,tmp,MINUS);
333 
334  r[rb[1]] -= b;
335  Double resid = sqrt(norm2(r,rb[1]));
336  Double resid_rel = resid/sqrt(norm2(b,rb[1]));
337  QDPIO::cout << "MdagM check: || r || = " << resid << " || r || / || b ||=" << resid_rel << std::endl;
338  ASSERT_LT( toDouble(resid_rel),1.0e-8);
339 
340 
341 
342 }
343 
344 #ifdef BUILD_QUDA
345 TEST_P(MdagMInvTestAsymm, CheckMdagMInvAsymmQUDAPredict)
346 {
347  std::istringstream inv_param_xml_stream(GetParam());
348  XMLReader xml_in(inv_param_xml_stream);
349 
350  GroupXML_t inv_param = readXMLGroup(xml_in, "//InvertParam", "invType");
351  Handle<MdagMSystemSolver<T>> MdagM_solver = S_asymm->invMdagM(state,inv_param);
352 
353  T b; gaussian(b,rb[1]);
354  T x; x[rb[1]] = zero;
355 
356  QUDA4DChronoPredictor chrono(5,DEFAULT);
357 
358  (*MdagM_solver)(x,b,chrono);
359 
360  T tmp, r;
361 
362  (*M_asymm)(tmp,x,PLUS);
363  (*M_asymm)(r,tmp,MINUS);
364 
365  r[rb[1]] -= b;
366  Double resid = sqrt(norm2(r,rb[1]));
367  Double resid_rel = resid/sqrt(norm2(b,rb[1]));
368  QDPIO::cout << "MdagM check: || r || = " << resid << " || r || / || b ||=" << resid_rel << std::endl;
369  ASSERT_LT( toDouble(resid_rel),2.0e-8);
370 
371 }
372 #endif
373 
374 INSTANTIATE_TEST_CASE_P(MdagMInvSysSolver,
375  MdagMInvTestSymm,
376  ::testing::Values(inv_param_syssolver_bicgstab_xml));
377 
378 #ifdef BUILD_QUDA
379 INSTANTIATE_TEST_CASE_P(MdagMInvQUDASolver,
380  MdagMInvTestSymm,
381  ::testing::Values(inv_param_quda_bicgstab_xml,
382  inv_param_quda_multigrid_xml));
383 #endif
384 
385 INSTANTIATE_TEST_CASE_P(MdagMInvSysSolver,
386  MdagMInvTestAsymm,
387  ::testing::Values(inv_param_syssolver_bicgstab_xml));
388 
389 #ifdef BUILD_QUDA
390 INSTANTIATE_TEST_CASE_P(MdagMInvQUDASolver,
391  MdagMInvTestAsymm,
392  ::testing::Values(inv_param_quda_bicgstab_asymm_xml,
393  inv_param_quda_multigrid_asymm_xml));
394 #endif
395 
396 class multiMdagMInvTestSymm : public SymmFixtureT<::testing::TestWithParam<std::string>>{};
397 class multiMdagMInvTestAsymm : public SymmFixtureT<::testing::TestWithParam<std::string>>{};
398 
399 TEST_P(multiMdagMInvTestSymm, checkMultiShift)
400 {
401  std::istringstream inv_param_xml_stream(GetParam());
402  XMLReader xml_in(inv_param_xml_stream);
403 
404  GroupXML_t inv_param = readXMLGroup(xml_in, "//InvertParam", "invType");
405 
406  Handle<MdagMMultiSystemSolver<T>> multiMdagM = S_symm->mInvMdagM(state,inv_param);
407 
408  T rhs;
409  gaussian(rhs,rb[1]);
410 
411  int n_shift=3;
412  multi1d<Real> shifts(n_shift);
413  shifts[0] = 0.0001;
414  shifts[1] = 0.01;
415  shifts[2] = 0.1;
416 
417  // Zero the initial guesses
418  multi1d<T> solns(n_shift);
419  for(int shift=0; shift < n_shift; ++shift) {
420  (solns[shift])[rb[1]]=zero;
421  }
422 
423  //operator() (multi1d< multi1d<T> >& psi, const multi1d<Real>& shifts, const multi1d<T>& chi)
424  (*multiMdagM)(solns,shifts,rhs);
425 
426  for(int shift = 0; shift < n_shift; ++shift) {
427  T r = zero;
428  T tmp = zero;
429  // r = M^\dag M solns[shift]
430  (*M_symm)(tmp,solns[shift],PLUS);
431  (*M_symm)(r, tmp, MINUS);
432 
433  // r = M^\dag M solns[shift] + shifts[shift] solns[shift]
434  // = (M^\dag M + shifts[shift]) solns[shift]
435  r[rb[1]] += shifts[shift]*solns[shift];
436 
437  // -residudum
438  r[rb[1]] -= rhs;
439 
440  Double resid_rel = sqrt(norm2(r,rb[1])/norm2(rhs,rb[1]));
441  QDPIO::cout << "shift="<<shift << " || r || / || b ||=" << resid_rel << std::endl;
442  ASSERT_LT( toDouble(resid_rel), 1.0e-8);
443  }
444 }
445 
446 TEST_P(multiMdagMInvTestAsymm, checkMultiShift)
447 {
448  std::istringstream inv_param_xml_stream(GetParam());
449  XMLReader xml_in(inv_param_xml_stream);
450 
451  GroupXML_t inv_param = readXMLGroup(xml_in, "//InvertParam", "invType");
452 
453  Handle<MdagMMultiSystemSolver<T>> multiMdagM = S_asymm->mInvMdagM(state,inv_param);
454 
455  T rhs;
456  gaussian(rhs,rb[1]);
457 
458  int n_shift=3;
459  multi1d<Real> shifts(n_shift);
460  shifts[0] = 0.0001;
461  shifts[1] = 0.01;
462  shifts[2] = 0.1;
463 
464  // Zero the initial guesses
465  multi1d<T> solns(n_shift);
466  for(int shift=0; shift < n_shift; ++shift) {
467  (solns[shift])[rb[1]]=zero;
468  }
469 
470  //operator() (multi1d< multi1d<T> >& psi, const multi1d<Real>& shifts, const multi1d<T>& chi)
471  (*multiMdagM)(solns,shifts,rhs);
472 
473  for(int shift = 0; shift < n_shift; ++shift) {
474  T r = zero;
475  T tmp = zero;
476  // r = M^\dag M solns[shift]
477  (*M_asymm)(tmp,solns[shift],PLUS);
478  (*M_asymm)(r, tmp, MINUS);
479 
480  // r = M^\dag M solns[shift] + shifts[shift] solns[shift]
481  // = (M^\dag M + shifts[shift]) solns[shift]
482  r[rb[1]] += shifts[shift]*solns[shift];
483 
484  // -residudum
485  r[rb[1]] -= rhs;
486 
487  Double resid_rel = sqrt(norm2(r,rb[1])/norm2(rhs,rb[1]));
488  QDPIO::cout << "shift="<<shift << " || r || / || b ||=" << resid_rel << std::endl;
489  ASSERT_LT( toDouble(resid_rel), 1.0e-8);
490  }
491 }
492 
493 INSTANTIATE_TEST_CASE_P(MultiShiftSysSolver,
494  multiMdagMInvTestAsymm,
495  ::testing::Values(inv_param_multi_cg_xml));
496 
497 #ifdef BUILD_QUDA
498 INSTANTIATE_TEST_CASE_P(MultiShiftQUDASolver,
499  multiMdagMInvTestAsymm,
500  ::testing::Values(inv_param_multi_cg_quda_asymm_xml));
501 #endif
502 
503 INSTANTIATE_TEST_CASE_P(MultiShiftSysSolver,
504  multiMdagMInvTestSymm,
505  ::testing::Values(inv_param_multi_cg_xml));
506 
507 #ifdef BUILD_QUDA
508 INSTANTIATE_TEST_CASE_P(MultiShiftQUDASolver,
509  multiMdagMInvTestSymm,
510  ::testing::Values(inv_param_multi_cg_quda_xml));
511 #endif
512 
513 
514 // Forces
515 TEST_F(SymmFixture, TestDeriv)
516 {
517  // M_symm = M_oo^{-1} M_asymm
518  // X^\dagger d[ M_symm ] Y
519  // = X^\dagger d[ M_oo^{-1} ] M_asymm Y
520  // +X^\dagger M_oo^{-1} d[ M_asymm ] Y
521  //
522  // = -X^\dagger M_oo^{-1} d[ M_oo ] M^{-1}_oo M_asymm Y
523  // +X^\dagger M_oo^{-1} d[ M_asymm ] Y
524  //
525  // = -Z^\dagger d[ M_oo ] W
526  // +Z^\dagger d[ M_asymm ] Y
527  //
528  // with W = M_{oo}^{-1} M_asymm Y
529  // and Z = M_oo^{-dagger} X
530 
531  P ds_symm;
532  P ds_tmp;
533  P rhs;
534 
535 
536  T X = zero;
537  T Y = zero;
538 
539  gaussian(X,rb[1]);
540  gaussian(Y,rb[1]);
541 
542  T tmp = zero;
543  T W = zero;
544  T Z = zero;
545 
546  // W = M^{-1}_oo M_asymm Y
547  (*M_asymm)(tmp,Y,PLUS);
548  (*M_symm).unprecOddOddInvLinOp(W,tmp,PLUS);
549 
550  (*M_symm).unprecOddOddInvLinOp(Z,X,MINUS);
551 
552  // The derivative of M_symm
553  (*M_symm).deriv(ds_symm,X,Y,PLUS);
554 
555  // rhs = Z^\dagger d[ M_asymm ] Y
556  (*M_asymm).deriv(rhs,Z,Y,PLUS);
557 
558  // rhs -= Z^\dagger d[ M_oo ] W
559  (*M_symm).derivUnprecOddOddLinOp(ds_tmp,Z,W,PLUS);
560  rhs -= ds_tmp;
561 
562  for(int mu=0; mu < Nd; ++mu) {
563  rhs[mu] -= ds_symm[mu];
564  Double norm_rhs = sqrt(norm2(rhs[mu]));
565  Double norm_rhs_per_number = norm_rhs/Double(3*3*2*Layout::vol());
566  QDPIO::cout << "mu=" << mu << " || rhs - ds_symm || = " << norm_rhs
567  << " || rhs - ds_symm || / number =" << norm_rhs_per_number << std::endl;
568 
569  ASSERT_LT(toDouble(norm_rhs_per_number), 1.0e-18 );
570  }
571 }
572 
573 TEST_F(SymmFixture, TestDerivDagger)
574 {
575  // M_symm^\dagger = M_asymm^\dagger M_oo^{-\dagger}
576 
577  // X^\dagger d[ M_symm ] Y
578  // = X^\dagger d[ M_asymm^\dagger ] M_oo^{-\dagger} Y
579  // +X^\dagger M_asymm^\dagger d[ M_oo^{-\dagger} ] Y
580  //
581  // = +X^\dagger d[M_asymm^\dagger ] M_oo^{-dagger} Y
582  // -X^\dagger M_asymm^\dagger ] M_oo^{-dagger} d[ M_oo^{\dagger} ] M_oo^{-dagger} Y
583  //
584  // = +X^\dagger d[ M_asymm^\dagger ] W
585  // -Z^\dagger d[ M_oo^\dagger ] W
586  //
587  // with W = M_{oo}^{-\dagger} Y
588  // and Z = M_oo^{-1} M_asymm X
589 
590  P ds_symm;
591  P ds_tmp;
592  P rhs;
593 
594 
595  T X = zero;
596  T Y = zero;
597 
598  gaussian(X,rb[1]);
599  gaussian(Y,rb[1]);
600 
601  T tmp = zero;
602  T W = zero;
603  T Z = zero;
604 
605  // W = M^{-dagger}_oo Y
606  (*M_symm).unprecOddOddInvLinOp(W,Y,MINUS);
607 
608  //Z = M_oo^{-1} M_asymm X
609  (*M_asymm)(tmp,X,PLUS);
610  (*M_symm).unprecOddOddInvLinOp(Z,tmp,PLUS);
611 
612  // The derivative of M_symm^dagger
613  (*M_symm).deriv(ds_symm,X,Y,MINUS);
614 
615  // rhs = X^\dagger d[ M_asymm ] W
616  (*M_asymm).deriv(rhs,X,W,MINUS);
617 
618  // rhs -= Z^\dagger d[ M_oo ] W
619  (*M_symm).derivUnprecOddOddLinOp(ds_tmp,Z,W,MINUS);
620  rhs -= ds_tmp;
621 
622  for(int mu=0; mu < Nd; ++mu) {
623  rhs[mu] -= ds_symm[mu];
624  Double norm_rhs = sqrt(norm2(rhs[mu]));
625  Double norm_rhs_per_number = norm_rhs/Double(3*3*2*Layout::vol());
626  QDPIO::cout << "mu=" << mu << " || rhs - ds_symm || = " << norm_rhs
627  << " || rhs - ds_symm || / number =" << norm_rhs_per_number << std::endl;
628 
629  ASSERT_LT(toDouble(norm_rhs_per_number), 1.0e-18 );
630  }
631 }
632 
633 TEST_F(SymmFixture,TestLogDetUnitGauge)
634 {
635  Q u_unit;
636  u_unit.resize(Nd);
637  for(int mu=0; mu < Nd; ++mu) {
638  u_unit[mu] = 1;
639  }
640 
641  Handle< FermState<T,P,Q> > unit_state = S_symm->createState(u_unit);
642  QDPIO::cout << "Unit Gauge Symm Linop creation" << std::endl;
643  Handle<LinOpSymm_T> lop = dynamic_cast<LinOpSymm_T*>(S_symm->linOp(unit_state));
644 
645  Double S_e = lop->logDetEvenEvenLinOp();
646  Double S_o = lop->logDetOddOddLinOp();
647  QDPIO::cout << "Unit gauge: log Det EvenEven = " << S_e << std::endl;
648  QDPIO::cout << "Unit gauge: log Det OddOdd = " << S_o << std::endl;
649 
650  Double absdiff = fabs(S_e-S_o);
651  QDPIO::cout << "Unit gauge: | S_e - S_o | = " << absdiff << std::endl;
652  ASSERT_LT( toDouble(absdiff), 1.0e-14);
653 }
654 
655 TEST_F(SymmFixture,TestLogDetShiftedGauge)
656 {
657  Q u_shifted;
658  u_shifted.resize(Nd);
659  for(int mu=0; mu < Nd; ++mu) {
660  u_shifted[mu] = shift(u[mu],FORWARD,3);
661  }
662 
663  Handle< FermState<T,P,Q> > shifted_state = S_symm->createState(u_shifted);
664 
665  QDPIO::cout << "Shifted Gauge Symm Linop creation" << std::endl;
666 
667  Handle<LinOpSymm_T> shifted_lop = dynamic_cast<LinOpSymm_T*>(S_symm->linOp(shifted_state));
668 
669  Double S_e_asymm = M_asymm->logDetEvenEvenLinOp();
670  Double S_e = M_symm->logDetEvenEvenLinOp();
671  Double S_o = M_symm->logDetOddOddLinOp();
672 
673  Double S_e_shifted = shifted_lop->logDetEvenEvenLinOp();
674  Double S_o_shifted = shifted_lop->logDetOddOddLinOp();
675 
676  QDPIO::cout << "Asymm op log Det EvenEven = " << S_e_asymm << std::endl;
677  QDPIO::cout << "Random gauge: log Det EvenEven = " << S_e << std::endl;
678  QDPIO::cout << "Random gauge: log Det OddOdd = " << S_o << std::endl;
679  QDPIO::cout << "Shifted gauge: log Det EvenEven = " << S_e_shifted << std::endl;
680  QDPIO::cout << "Shifted gauge: log Det OddOdd = " << S_o_shifted << std::endl;
681 
682  Double diff_e_symm_asymm = fabs(S_e_asymm - S_e);
683  Double diff_eo = fabs(S_e - S_o_shifted)/Double(rb[0].numSiteTable());
684  Double diff_oe = fabs(S_o - S_e_shifted)/Double(rb[1].numSiteTable());
685  QDPIO::cout << "| logDet_e_asymm - logdet_e_symm | = " << diff_e_symm_asymm << std::endl;
686  QDPIO::cout << "| logDet_e - logDet_o_shifted |/site = " << diff_eo << std::endl;
687  QDPIO::cout << "| logDet_o - logDet_e_shifted |/site = " << diff_oe << std::endl;
688 
689  ASSERT_LT( toDouble(diff_e_symm_asymm), 1.0e-14);
690  ASSERT_LT( toDouble(diff_eo), 1.0e-13);
691  ASSERT_LT( toDouble(diff_oe), 1.0e-13);
692 
693 }
694 
695 TEST_F(SymmFixture,TestTwist)
696 {
697 
698  Real twist=Real(0.05);
699 
700 
701 
702 
703  LatticeFermion source;
704  gaussian(source,rb[1]);
705  LatticeFermion t1, t2;
706 
707  {
708  (*M_tw)(t1,source,PLUS);
709  (*M_symm)(t2,source, PLUS);
710  t2[ rb[1] ] += twist*(GammaConst<Ns,Ns*Ns-1>()*timesI(source));
711 
712  t2[ rb[1] ] -= t1;
713  Double normdiff = sqrt(norm2(t2,rb[1]));
714  QDPIO::cout << "PLUS : || M(mu) - ( Mdag + i gamma_5 mu ) || = "
715  << normdiff << std::endl;
716 
717  ASSERT_LT( toDouble(normdiff), 1.0e-14);
718  }
719 
720  {
721  (*M_tw)(t1,source,MINUS);
722  (*M_symm)(t2,source, MINUS);
723  t2[ rb[1] ] -= twist*(GammaConst<Ns,Ns*Ns-1>()*timesI(source));
724 
725  t2[ rb[1] ] -= t1;
726  Double normdiff = sqrt(norm2(t2,rb[1]));
727  QDPIO::cout << "MINUS : || M^dag(mu) - ( Mdag - igamma5 mu ) || = " << normdiff << std::endl;
728 
729  ASSERT_LT( toDouble(normdiff), 1.0e-14);
730  }
731 
732  {
733  LatticeFermion mdagm,t3;
734  (*M_tw)(t1,source, PLUS);
735  (*M_tw)(t2,t1,MINUS);
736 
737  // M^\dag M
738  (*M_symm)(t1,source,PLUS);
739  (*M_symm)(mdagm,t1,MINUS);
740 
741  // + mu^2 source
742  mdagm[rb[1]] += (twist*twist)*source;
743 
744  // +i mu M^\dag gamma_5
745  t1[rb[1]] = (GammaConst<Ns,Ns*Ns-1>()*timesI(source));
746  (*M_symm)(t3,t1,MINUS);
747  mdagm[rb[1]] += twist*t3;
748 
749  // -i mu gamma_5 M soure
750  (*M_symm)(t1,source,PLUS);
751  t3[rb[1]] = (GammaConst<Ns,Ns*Ns-1>()*timesI(t1));
752  mdagm[rb[1]] -= twist*t3;
753 
754 
755  mdagm[rb[1]] -= t2;
756  Double normdiff = sqrt(norm2(mdagm,rb[1]));
757  QDPIO::cout << "MDAGM : || M^dag(mu)M(mu) - ( Mdag M + mu^2 + imu[ M^dag g_5 - g_5 M ]) || = " << normdiff << std::endl;
758  ASSERT_LT( toDouble(normdiff), 1.0e-13);
759  }
760 }
761 
762 class TrLogForceFixture : public SymmFixtureT<::testing::TestWithParam<enum PlusMinus>>{};
763 TEST_P(TrLogForceFixture,TestShiftedGaugeTrLnForce)
764 {
765  Q u_shifted;
766  u_shifted.resize(Nd);
767  for(int mu=0; mu < Nd; ++mu) {
768  u_shifted[mu] = shift(u[mu],FORWARD,3);
769  }
770 
771  Handle< FermState<T,P,Q> > periodic_state = S_symm_periodic->createState(u);
772  Handle< FermState<T,P,Q> > periodic_shifted_state = S_symm_periodic->createState(u_shifted);
773 
774  QDPIO::cout << "Shifted Gauge Symm Linop creation" << std::endl;
775 
776  Handle<LinOpAsymm_T> asymm_periodic_op = dynamic_cast<LinOpAsymm_T*>(S_asymm_periodic->linOp(periodic_state));
777  Handle<LinOpSymm_T> symm_periodic_op = dynamic_cast<LinOpSymm_T*>(S_symm_periodic->linOp(periodic_state));
778  Handle<LinOpSymm_T> shifted_periodic_op = dynamic_cast<LinOpSymm_T*>(S_symm_periodic->linOp(periodic_shifted_state));
779 
780  //! Get the force from the EvenEven Trace Log
781  // virtual void derivLogDetEvenEvenLinOp(P& ds_u, enum PlusMinus isign) const = 0;
782  P ds_asymm;
783  asymm_periodic_op->derivLogDetEvenEvenLinOp(ds_asymm,GetParam());
784 
785  P ds_symm;
786  symm_periodic_op->derivLogDetEvenEvenLinOp(ds_symm,GetParam());
787 
788  P ds_symm_shifted;
789  shifted_periodic_op->derivLogDetOddOddLinOp(ds_symm_shifted, GetParam());
790 
791  P ds_unshifted; ds_unshifted.resize(Nd);
792  P ds_wrong_shifted; ds_wrong_shifted.resize(Nd);
793  for(int mu=0; mu < Nd; ++mu ) {
794  ds_unshifted[mu] = shift(ds_symm_shifted[mu], BACKWARD,3);
795  ds_wrong_shifted[mu] = shift(ds_symm_shifted[mu], FORWARD, 3);
796  }
797  for(int mu=0; mu < Nd; ++mu ) {
798  Double mu_contrib_asymm = norm2(ds_asymm[mu]);
799 
800  Double mu_contrib_symm = norm2(ds_symm[mu]);
801 
802  Double mu_contrib_shifted = norm2(ds_symm_shifted[mu]);
803 
804  QDPIO::cout << "mu=" << mu << " asymm force norm="<< mu_contrib_asymm << std::endl;
805  QDPIO::cout << "mu=" << mu << " symm force norm ="<< mu_contrib_symm << std::endl;
806  QDPIO::cout << "mu=" << mu << " shifted force norm="<< mu_contrib_shifted << std::endl;
807 
808 
809  Double diff_symm_asymm_ee = fabs(mu_contrib_asymm - mu_contrib_symm );
810  Double diff_symm_ee_oo = fabs(mu_contrib_symm - mu_contrib_shifted );
811 
812  QDPIO::cout << "mu=" << mu << " | F_asymm_ee - F_symm_ee | = " << diff_symm_asymm_ee << std::endl;
813  QDPIO::cout << "mu=" << mu << " | F_ee - F_shifted_oo | = " << diff_symm_ee_oo << std::endl;
814 
815  ASSERT_LT( toDouble(diff_symm_asymm_ee), 1.0e-15);
816  ASSERT_LT( toDouble(diff_symm_ee_oo), 1.0e-15);
817 
818  ds_unshifted[mu] -= ds_symm[mu];
819  ds_wrong_shifted[mu] -= ds_symm[mu];
820  Double diffnorm_unshifted = sqrt(norm2(ds_unshifted[mu]));
821  Double diffnorm_per_site = diffnorm_unshifted/Layout::vol();
822  QDPIO::cout << "mu=" << mu << " || F_mu - shift(F_shifted, BACKWARD, mu) || =" << diffnorm_unshifted << std::endl;
823  QDPIO::cout << "mu=" << mu << " || F_mu - shift(F_shifted, BACKWARD, mu) ||/site =" << diffnorm_per_site << std::endl;
824 
825  ASSERT_LT( toDouble(diffnorm_per_site), 1.0e-14 );
826 
827  Double diffnorm_wrong_shifted = sqrt(norm2(ds_wrong_shifted[mu]));
828  QDPIO::cout << "mu=" << mu << " DELIBERATELY WRONG: || F_mu - shift(F_shifted, FORWARD, mu) || =" << diffnorm_wrong_shifted << std::endl;
829  ASSERT_GT( toDouble(diffnorm_wrong_shifted), 0.5);
830  QDPIO::cout << std::endl;
831  }
832 
833 }
834 
835 INSTANTIATE_TEST_CASE_P(TrLogForces,
836  TrLogForceFixture,
837  ::testing::Values(PLUS,MINUS));
838 
839 
840 TEST_F(SymmFixture, CheckDerivMultipole)
841 {
842  int n_poles = 5;
843  multi1d<T> X(n_poles);
844  multi1d<T> Y(n_poles);
845 
846 
847  for(int i=0; i < n_poles; ++i) {
848  gaussian(X[i],rb[1]);
849  gaussian(Y[i],rb[1]);
850  }
851 
852  multi1d<LatticeColorMatrix> ds_u(Nd);
853  multi1d<LatticeColorMatrix> ds_tmp(Nd);
854  multi1d<LatticeColorMatrix> ds_mp(Nd);
855 
856  for(int dag=0; dag < 2; ++dag ) {
857  QDPIO::cout << "Dagger = " << dag << std::endl;
858 
859  enum PlusMinus isign = (dag == 0 ) ? PLUS : MINUS;
860  // Isign = PLUS
861  // Zero sum
862  for(int mu=0; mu < Nd; ++mu) {
863  ds_u[mu] = zero;
864  }
865 
866  // Accumulate forces
867  StopWatch swatch;
868  swatch.start();
869  for(int i=0; i < n_poles; ++i) {
870  M_symm->deriv(ds_tmp, X[i],Y[i], isign);
871  ds_u += ds_tmp;
872  }
873  swatch.stop();
874  QDPIO::cout << "Individual Poles took: " << swatch.getTimeInSeconds() << " sec." << std::endl;
875 
876  // do Multipole version
877  swatch.reset(); swatch.start();
878  M_symm->derivMultipole(ds_mp, X,Y,isign);
879  swatch.stop();
880  QDPIO::cout << "Multipole took: " << swatch.getTimeInSeconds() << " sec." << std::endl;
881 
882  for(int mu=0; mu < Nd; ++mu) {
883  ds_mp[mu] -= ds_u[mu];
884  Double normdiff = sqrt(norm2(ds_mp[mu]));
885  Double normdiff_per_link = normdiff/static_cast<double>(4*Layout::vol());
886  QDPIO::cout << "mu="<< mu << " normdiff = " << normdiff
887  << " normdiff per link = " << normdiff_per_link << std::endl;
888  ASSERT_LT( toDouble(normdiff_per_link), 1.0e-17);
889  }
890  }
891 }
chromabase.h
Primary include file for CHROMA library code.
Chroma::EvenOddPrecCloverLinOp
Even-odd preconditioned Clover-Dirac operator.
Definition: eoprec_clover_linop_w.h:31
Chroma::EvenOddPrecWilsonTypeFermAct
Even-odd preconditioned Wilson-like fermion actions including derivatives.
Definition: eoprec_wilstype_fermact_w.h:23
Chroma::QUDA4DChronoPredictor
Zero initial guess predictor.
Definition: quda_predictor.h:32
Chroma::SingletonHolder::Instance
static T & Instance()
Definition: singleton.h:432
Chroma::SymEvenOddPrecCloverLinOp
Symmetric even-odd preconditioned Clover-Dirac operator.
Definition: seoprec_clover_linop_w.h:31
Chroma::SymEvenOddPrecWilsonTypeFermAct
Symmetric even-odd preconditioned Wilson-like fermion actions including derivatives.
Definition: seoprec_wilstype_fermact_w.h:22
SymmFixtureT::S_symm_periodic
Handle< S_symm_T > S_symm_periodic
Definition: symm_prec_tests.cc:94
SymmFixtureT::M_symm
Handle< LinOpSymm_T > M_symm
Definition: symm_prec_tests.cc:100
SymmFixtureT::state
Handle< FermState< T, P, Q > > state
Definition: symm_prec_tests.cc:97
SymmFixtureT::S_asymm
Handle< S_asymm_T > S_asymm
Definition: symm_prec_tests.cc:93
SymmFixtureT::S_asymm_periodic
Handle< S_asymm_T > S_asymm_periodic
Definition: symm_prec_tests.cc:95
SymmFixtureT::M_tw
Handle< LinOpSymm_T > M_tw
Definition: symm_prec_tests.cc:101
SymmFixtureT::Q
multi1d< LatticeColorMatrix > Q
Definition: symm_prec_tests.cc:31
SymmFixtureT::M_asymm
Handle< LinOpAsymm_T > M_asymm
Definition: symm_prec_tests.cc:99
SymmFixtureT::S_symm
Handle< S_symm_T > S_symm
Definition: symm_prec_tests.cc:92
SymmFixtureT::T
LatticeFermion T
Definition: symm_prec_tests.cc:30
SymmFixtureT::S_symm_twisted
Handle< S_symm_T > S_symm_twisted
Definition: symm_prec_tests.cc:96
SymmFixtureT::P
multi1d< LatticeColorMatrix > P
Definition: symm_prec_tests.cc:32
mu
int mu
Definition: cool.cc:24
eoprec_clover_linop_w.h
Even-odd preconditioned Clover fermion linear operator.
eoprec_linop.h
Base class for even-odd preconditioned 4D and 5D Linop.
eoprec_wilstype_fermact_w.h
Even-odd preconditioned Wilson-like fermion actions.
fermact_factory_w.h
Fermion action factories.
Chroma::readXMLGroup
GroupXML_t readXMLGroup(XMLReader &xml_in, const std::string &path, const std::string &type_name)
Read group and return as a std::string.
Definition: xml_group_reader.cc:11
ASSERT_GT
#define ASSERT_GT(val1, val2)
Definition: gtest.h:19801
ASSERT_LT
#define ASSERT_LT(val1, val2)
Definition: gtest.h:19793
handle.h
Class for counted reference semantics.
x
int x
Definition: meslate.cc:34
Nd
Nd
Definition: meslate.cc:74
Chroma
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
Chroma::gaussian
gaussian(aux)
Chroma::u
static multi1d< LatticeColorMatrix > u
Definition: syssolver_linop_qop_mg_w.cc:39
Chroma::Q
LinOpSysSolverMGProtoClover::Q Q
Definition: syssolver_linop_clover_mg_proto.cc:64
Chroma::tmp
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
Chroma::T
LinOpSysSolverMGProtoClover::T T
Definition: syssolver_linop_clover_mg_proto.cc:63
Chroma::i
int i
Definition: pbg5p_w.cc:55
Chroma::MINUS
@ MINUS
Definition: chromabase.h:45
Chroma::PLUS
@ PLUS
Definition: chromabase.h:45
Chroma::reunit
void reunit(LatticeColorMatrixF3 &xa)
Definition: reunit.cc:467
Chroma::b
Complex b
Definition: invbicg.cc:96
Chroma::state
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > state
Definition: pbg5p_w.cc:28
Chroma::zero
Double zero
Definition: invbicg.cc:106
SymmPrecTesting::fermact_xml_symm
std::string fermact_xml_symm
Definition: symm_prec_xml.h:41
SymmPrecTesting::inv_param_multi_cg_xml
std::string inv_param_multi_cg_xml
Definition: symm_prec_xml.h:367
SymmPrecTesting::fermact_xml_asymm_periodic
std::string fermact_xml_asymm_periodic
Definition: symm_prec_xml.h:92
SymmPrecTesting::fermact_xml_asymm
std::string fermact_xml_asymm
Definition: symm_prec_xml.h:16
SymmPrecTesting::inv_param_syssolver_bicgstab_xml
std::string inv_param_syssolver_bicgstab_xml
Definition: symm_prec_xml.h:144
SymmPrecTesting::fermact_xml_symm_periodic
std::string fermact_xml_symm_periodic
Definition: symm_prec_xml.h:118
SymmPrecTesting::fermact_xml_symm_twisted
std::string fermact_xml_symm_twisted
Definition: symm_prec_xml.h:65
testing::internal::Double
FloatingPoint< double > Double
Definition: gtest.h:7351
testing::Values
internal::ValueArray1< T1 > Values(T1 v1)
Definition: gtest.h:15914
FORWARD
#define FORWARD
Definition: primitives.h:82
BACKWARD
#define BACKWARD
Definition: primitives.h:83
quda_predictor.h
Pick channel for QUDA Predictor.
reunit.h
Reunitarize in place a color matrix to SU(N)
seoprec_clover_linop_w.h
Symmetric even-odd preconditioned Clover fermion linear operator.
seoprec_linop.h
Base class for symmetric even-odd preconditioned 4D and 5D Linop.
seoprec_wilstype_fermact_w.h
Symmetric even-odd preconditioned Wilson-like fermion actions.
Chroma::GroupXML_t
Hold group xml and type id.
Definition: xml_group_reader.h:17
TEST_P
TEST_P(QPropTest, CheckQprop)
Definition: symm_prec_tests.cc:210
TEST_F
TEST_F(SymmFixture, CheckOp)
Definition: symm_prec_tests.cc:107
INSTANTIATE_TEST_CASE_P
INSTANTIATE_TEST_CASE_P(PropSyssolver, QPropTest, ::testing::Values(inv_param_syssolver_bicgstab_xml))
symm_prec_xml.h
P
multi1d< LatticeColorMatrix > P
Definition: t_clover.cc:13
xml_group_reader.h
Read an XML group as a std::string.
Generated by doxygen 1.9.1