CHROMA
t_temp_prec.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Test 4d fermion actions
3  */
4 
5 #include <iostream>
6 #include <cstdio>
7 
8 #include "chroma.h"
9 
10 
11 using namespace Chroma;
12 
13 #include "actions/ferm/linop/eo3dprec_s_cprec_t_wilson_linop_w.h"
14 #include "actions/ferm/linop/eo3dprec_s_cprec_t_clover_linop_w.h"
15 //! To insure linking of code, place the registered code flags here
16 /*! This is the bit of code that dictates what fermacts are in use */
17 bool linkage_hack()
18 {
19  bool foo = true;
20  // All actions
21  foo &= WilsonTypeFermActsEnv::registerAll();
22  return foo;
23 }
24 
25 
26 
27 int main(int argc, char **argv)
28 {
29  // This test has trouble with Nc=2, so make sure we're using 3 colors
30 #if QDP_NC == 3
31  // Put the machine into a known state
32  Chroma::initialize(&argc, &argv);
33 
34  QDPIO::cout << "linkage=" << linkage_hack() << std::endl;
35 
36 
37  multi1d<int> nrow(Nd);
38  for(int i=0; i < Nd-1; i++) {
39  nrow[i] = 4;
40  }
41  nrow[Nd-1] = 16;
42 
43  // Specify lattice size, shape, etc.
44  Layout::setLattSize(nrow);
45  Layout::create();
46 
47  QDPIO::cout << "t_temp_prec" << std::endl;
48 
49  // Start up a weak field
50  struct Cfg_t config = { CFG_TYPE_WEAK_FIELD, "dummy" };
51  multi1d<LatticeColorMatrix> u(Nd);
52  XMLReader gauge_file_xml, gauge_xml;
53  gaugeStartup(gauge_file_xml, gauge_xml, u, config);
54  // Check if the gauge field configuration is unitarized
55  unitarityCheck(u);
56 
57  // Instantiate XML writer for XMLDAT
58  XMLFileWriter& xml_out = Chroma::getXMLOutputInstance();
59  push(xml_out, "t_temp_prec");
60  proginfo(xml_out); // Print out basic program info
61 
62 
63  // Calculate some gauge invariant observables just for info.
64  MesPlq(xml_out, "Observables", u);
65  xml_out.flush();
66  pop(xml_out);
67 
68 #if 0
69  // Now need to create a SimpleFermState< LatticeFermion,
70  // multi1d<LatticeColorMatrix>,
71  // multi1d<LatticeColorMatrix> >
72  multi1d<int> boundary(Nd);
73  boundary[0] = 1;
74  boundary[1] = 1;
75  boundary[2] = 1;
76  boundary[3] = -1;
77 
78 
79 #endif
80 
81  XMLReader xml_in("./t_temp_prec.ini.xml");
82 
83  std::string fermact, prec_fermact, ilu_fermact;
84  read(xml_in, "/tempPrec/UnprecFermAct/FermAct", fermact);
85  read(xml_in, "/tempPrec/PrecFermAct/FermAct", prec_fermact);
86  read(xml_in, "/tempPrec/ILUPrecFermAct/FermAct", ilu_fermact);
87  // Typedefs to save typing
88  typedef LatticeFermion T;
89  typedef multi1d<LatticeColorMatrix> P;
90  typedef multi1d<LatticeColorMatrix> Q;
91 
92 
93  // Create unprec and Prec actions
94  Handle< FermionAction<T,P,Q> >
95  S_unprec( TheFermionActionFactory::Instance().createObject(fermact,
96  xml_in,
97  "/tempPrec/UnprecFermAct") );
98 
99  UnprecWilsonFermAct& S_un_wils = dynamic_cast< UnprecWilsonFermAct& >(*S_unprec);
100 
101  Handle< FermionAction<T,P,Q> >
102  S_prec( TheFermionActionFactory::Instance().createObject(prec_fermact,
103  xml_in,
104  "/tempPrec/PrecFermAct") );
105 
106  UnprecSpaceCentralPrecTimeWilsonFermAct& S_unprec_s_cprec_t =
107  dynamic_cast< UnprecSpaceCentralPrecTimeWilsonFermAct& >(*S_prec);
108 
109 
110  Handle< FermionAction<T,P,Q> >
111  S_iluprec_handle( TheFermionActionFactory::Instance().createObject(ilu_fermact,
112  xml_in,
113  "/tempPrec/ILUPrecFermAct") );
114 
115  ILUPrecSpaceCentralPrecTimeWilsonTypeFermAct<T,P,Q>& S_iluprec=
116  dynamic_cast< ILUPrecSpaceCentralPrecTimeWilsonTypeFermAct<T,P,Q>& >(*S_iluprec_handle);
117 
118 
119 
120  // Now create a FermState
121  Handle< FermState<T, P, Q> > fs( S_un_wils.createState(u) );
122 
123 
124  Handle< UnprecSpaceCentralPrecTimeLinearOperator<T,P,Q> > D_temp_prec_handle( S_unprec_s_cprec_t.linOp( fs ) );
125  UnprecSCprecTWilsonLinOp& D_temp_prec = dynamic_cast<UnprecSCprecTWilsonLinOp&>(*D_temp_prec_handle);
126 
127  Handle< UnprecLinearOperator<T,P,Q> > D_w_handle( S_un_wils.linOp(fs) );
128  UnprecWilsonLinOp& D_w = dynamic_cast< UnprecWilsonLinOp& >(*D_w_handle);
129 
130  Handle< ILUPrecSpaceCentralPrecTimeLinearOperator<T,P,Q> > D_tprec2_handle( S_iluprec.linOp(fs) );
131  ILUPrecSCprecTWilsonLinOp& D_temp_prec2 = dynamic_cast< ILUPrecSCprecTWilsonLinOp& >(*D_tprec2_handle);
132 
133  LatticeFermion chi;
134  LatticeFermion psi1, psi2, tmp1, tmp2;
135  LatticeFermion diff;
136 
137  psi1 = zero;
138  psi2 = zero;
139 
140  gaussian(chi);
141  D_temp_prec.getFermBC().modifyF(chi);
142 
143  // gamma_5 hermiticity: psi1 = \gamma_5 C_R^{-1} \gamma_5 \chi
144  tmp1 = Gamma(15)*chi;
145  D_temp_prec.invCRightLinOp(tmp2, tmp1, PLUS);
146  psi1 = Gamma(15)*tmp2;
147 
148  // psi2 = ( C_L^{-1} )^\dagger \chi
149  D_temp_prec.invCLeftLinOp(psi2, chi, MINUS);
150  diff = psi2 - psi1;
151  QDPIO::cout << " Gamma5_1 = " << sqrt( norm2(diff) / norm2(psi1) ) << std::endl;
152 
153 
154  // Other way
155  // psi1 = \gamma_5 C_L^{-1} \gamma_5 \chi
156  tmp1 = Gamma(15)*chi;
157  D_temp_prec.invCLeftLinOp(tmp2, tmp1, PLUS);
158  psi1 = Gamma(15)*tmp2;
159 
160  // psi2 = ( C_R^{-1} )^\dagger \chi
161  D_temp_prec.invCRightLinOp(psi2, chi, MINUS);
162  diff = psi2 - psi1;
163  QDPIO::cout << " Gamma5_2 = " << sqrt( norm2(diff) / norm2(psi1) ) << std::endl;
164 
165  // Check D_w = CL^{-1} CR^{-1} + D_s
166  // D_w
167  D_w(psi1, chi, PLUS);
168 
169  // CL^{-1} CR^{-1} + D_s)
170  D_temp_prec.invCRightLinOp(tmp1, chi, PLUS);
171  D_temp_prec.invCLeftLinOp(tmp2, tmp1, PLUS);
172  D_temp_prec.spaceLinOp(psi2, chi, PLUS);
173  psi2+=tmp2;
174 
175  diff = psi2 - psi1;
176  QDPIO::cout << " D_1 + = " << sqrt( norm2(diff)/norm2(psi1) ) << std::endl;
177 
178 
179  // Check D_w^\dag = CR^{-dagger} CL^{-dagger} + D_s^\dag (MINUS)
180  // D_w
181  D_w(psi1, chi, MINUS);
182 
183  // CL^{-1} CR^{-1} + D_s (MINUS)
184  D_temp_prec.invCLeftLinOp(tmp1, chi, MINUS);
185  D_temp_prec.invCRightLinOp(tmp2, tmp1, MINUS);
186  D_temp_prec.spaceLinOp(psi2, chi, MINUS);
187  psi2+=tmp2;
188 
189  diff = psi2 - psi1;
190  QDPIO::cout << " D_1 - = " << sqrt( norm2(diff)/norm2(psi1) ) << std::endl;
191 
192 #if 0
193  multi1d<LatticeColorMatrix> u_2(1);
194  multi1d<LatticeColorMatrix> u_inv(1);
195 
196  // Hotstart
197  HotSt(u_2);
198  PScalar<PColorMatrix< RComplex<REAL>, 3> > munit;
199  munit = Real(1).elem();
200 
201 
202  // This used to test the invert 3by3 - I made that private now.
203  // Tested through context
204 
205  munit = Real(1).elem(); // Invert
206  for(int site=all.start(); site <= all.end(); site++) {
207  u_2[0].elem(site) += munit;
208  D_temp_prec.invert3by3( u_inv[0].elem(site), u_2[0].elem(site) );
209  }
210 
211  LatticeColorMatrix identityP = u_2[0]*u_inv[0];
212  LatticeColorMatrix ident;
213  for(int site=all.start(); site <= all.end(); site++) {
214  identityP.elem(site) -= munit;
215  ident.elem(site) = munit;
216  }
217 
218  QDPIO::cout << "|| I - u_2*inv ||/ || 1 || = " << sqrt( norm2( identityP ) / norm2(ident) ) << std::endl;
219 #endif
220 
221 
222 #if 0
223  // Used to test the actual T0 and InvT0 ops - these are private now.
224  LatticeHalfFermion chi_half;
225  LatticeHalfFermion psi_half;
226  LatticeHalfFermion chi2_half;
227  LatticeHalfFermion diff_half;
228 
229  gaussian(chi_half);
230  D_temp_prec.TOp(psi_half, chi_half, PLUS);
231  D_temp_prec.invTOp(chi2_half, psi_half, PLUS);
232 
233  diff_half = chi2_half - chi_half;
234  QDPIO::cout << " || diff || / || chi_half || = " << sqrt(norm2(diff_half) / norm2(chi_half)) << std::endl;
235 
236 
237  D_temp_prec.TOp(psi_half, chi_half, MINUS);
238  D_temp_prec.invTOp(chi2_half, psi_half, MINUS);
239 
240  diff_half = chi2_half - chi_half;
241  QDPIO::cout << " || diff || / || chi_half || = " << sqrt(norm2(diff_half) / norm2(chi_half)) << std::endl;
242 #endif
243 
244  // Test LeftInverse is inverse of Left (Both orderings, PLUS and MINUS)
245  gaussian(chi);
246  D_temp_prec.getFermBC().modifyF(chi);
247 
248  D_temp_prec.cLeftLinOp(psi1, chi, PLUS);
249  D_temp_prec.invCLeftLinOp(tmp1, psi1, PLUS);
250 
251  diff = tmp1 - chi;
252  QDPIO::cout << " CL CL_INV + = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
253 
254 
255  D_temp_prec.cLeftLinOp(psi1, chi, MINUS);
256  D_temp_prec.invCLeftLinOp(tmp1, psi1, MINUS);
257 
258  diff = tmp1 - chi;
259  QDPIO::cout << " CL CL_INV - = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
260 
261  D_temp_prec.invCLeftLinOp(psi1, chi, PLUS);
262  D_temp_prec.cLeftLinOp(tmp1, psi1, PLUS);
263 
264  diff = tmp1 - chi;
265  QDPIO::cout << " CL_INV CL + = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
266 
267  D_temp_prec.invCLeftLinOp(psi1, chi, MINUS);
268  D_temp_prec.cLeftLinOp(tmp1, psi1, MINUS);
269 
270  diff = tmp1 - chi;
271  QDPIO::cout << " CL_INV CL - = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
272 
273 
274  // Test C_R is inverse of C_R_Inverse both orders, PLUS & Minus
275  gaussian(chi);
276  D_temp_prec.getFermBC().modifyF(chi);
277  D_temp_prec.cRightLinOp(psi1, chi, PLUS);
278  D_temp_prec.invCRightLinOp(tmp1, psi1, PLUS);
279 
280  diff = tmp1 - chi;
281  QDPIO::cout << " CR CR_INV + = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
282 
283  D_temp_prec.cRightLinOp(psi1, chi, MINUS);
284  D_temp_prec.invCRightLinOp(tmp1, psi1, MINUS);
285 
286  diff = tmp1 - chi;
287  QDPIO::cout << " CR CR_INV - = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
288 
289  D_temp_prec.invCRightLinOp(psi1, chi, PLUS);
290  D_temp_prec.cRightLinOp(tmp1, psi1, PLUS);
291 
292  diff = tmp1 - chi;
293  QDPIO::cout << " CR_INV CR + = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
294 
295  D_temp_prec.invCRightLinOp(psi1, chi, MINUS);
296  D_temp_prec.cRightLinOp(tmp1, psi1, MINUS);
297 
298  diff = tmp1 - chi;
299  QDPIO::cout << " CR_INV CR - = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
300 
301 
302  // Now test against the normal D_op = C_L^{-1}( 1 + C^L D_s C^R ) C_R^{-1}
303  D_w(psi1, chi, PLUS);
304  D_temp_prec.unprecLinOp(psi2, chi, PLUS);
305  diff = psi2 - psi1;
306  QDPIO::cout << " D_2 + = " << sqrt( norm2(diff) / norm2(psi1) ) << std::endl;
307 
308  D_w(psi1, chi, MINUS);
309  D_temp_prec.unprecLinOp(psi2, chi, MINUS);
310  diff = psi2 - psi1;
311  QDPIO::cout << " D_2 - = " << sqrt( norm2(diff) / norm2(psi1) ) << std::endl;
312 
313 
314  // ==============================================================================================
315  // ILU Prec Op
316  // ===========================================================================================
317 
318 
319  QDPIO::cout << "ILU Prec Op: " << std::endl;
320  // ILUPrecSCprecTWilsonLinOp D_temp_prec2(fs, Mass, aniso );
321 
322  tmp1 = Gamma(15)*chi;
323  D_temp_prec2.invCRightLinOp(tmp2, tmp1, PLUS);
324  psi1 = Gamma(15)*tmp2;
325 
326  // psi2 = ( C_L^{-1} )^\dagger \chi
327  D_temp_prec2.invCLeftLinOp(psi2, chi, MINUS);
328  diff = psi2 - psi1;
329  QDPIO::cout << " Gamma5_1 = " << sqrt( norm2(diff) / norm2(psi1) ) << std::endl;
330 
331 
332  // Other way
333  // psi1 = \gamma_5 C_L^{-1} \gamma_5 \chi
334  tmp1 = Gamma(15)*chi;
335  D_temp_prec2.invCLeftLinOp(tmp2, tmp1, PLUS);
336  psi1 = Gamma(15)*tmp2;
337 
338  // psi2 = ( C_R^{-1} )^\dagger \chi
339  D_temp_prec2.invCRightLinOp(psi2, chi, MINUS);
340  diff = psi2 - psi1;
341  QDPIO::cout << " Gamma5_2 = " << sqrt( norm2(diff) / norm2(psi1) ) << std::endl;
342 
343  // Test LeftInverse is inverse of Left (Both orderings, PLUS and MINUS)
344  gaussian(chi);
345  D_temp_prec2.getFermBC().modifyF(chi);
346  D_temp_prec2.cLeftLinOp(psi1, chi, PLUS);
347  D_temp_prec2.invCLeftLinOp(tmp1, psi1, PLUS);
348 
349  diff = tmp1 - chi;
350  QDPIO::cout << " CL CL_INV + = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
351 
352 
353  D_temp_prec2.cLeftLinOp(psi1, chi, MINUS);
354  D_temp_prec2.invCLeftLinOp(tmp1, psi1, MINUS);
355 
356  diff = tmp1 - chi;
357  QDPIO::cout << " CL CL_INV - = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
358 
359  D_temp_prec2.invCLeftLinOp(psi1, chi, PLUS);
360  D_temp_prec2.cLeftLinOp(tmp1, psi1, PLUS);
361 
362  diff = tmp1 - chi;
363  QDPIO::cout << " CL_INV CL + = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
364 
365  D_temp_prec2.invCLeftLinOp(psi1, chi, MINUS);
366  D_temp_prec2.cLeftLinOp(tmp1, psi1, MINUS);
367 
368  diff = tmp1 - chi;
369  QDPIO::cout << " CL_INV CL - = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
370 
371 
372  // Test C_R is inverse of C_R_Inverse both orders, PLUS & Minus
373  gaussian(chi);
374  D_temp_prec2.getFermBC().modifyF(chi);
375  D_temp_prec2.cRightLinOp(psi1, chi, PLUS);
376  D_temp_prec2.invCRightLinOp(tmp1, psi1, PLUS);
377 
378  diff = tmp1 - chi;
379  QDPIO::cout << " CR CR_INV + = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
380 
381  D_temp_prec2.cRightLinOp(psi1, chi, MINUS);
382  D_temp_prec2.invCRightLinOp(tmp1, psi1, MINUS);
383 
384  diff = tmp1 - chi;
385  QDPIO::cout << " CR CR_INV - = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
386 
387  D_temp_prec2.invCRightLinOp(psi1, chi, PLUS);
388  D_temp_prec2.cRightLinOp(tmp1, psi1, PLUS);
389 
390  diff = tmp1 - chi;
391  QDPIO::cout << " CR_INV CR + = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
392 
393  D_temp_prec2.invCRightLinOp(psi1, chi, MINUS);
394  D_temp_prec2.cRightLinOp(tmp1, psi1, MINUS);
395 
396  diff = tmp1 - chi;
397  QDPIO::cout << " CR_INV CR - = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
398 
399 
400  // Now test against the normal D_op = C_L^{-1} unprecOp() C_R^{-1}
401  D_w(psi1, chi, PLUS);
402  D_temp_prec2.unprecLinOp(psi2, chi, PLUS);
403  diff = psi2 - psi1;
404  QDPIO::cout << " D_2 + = " << sqrt( norm2(diff) / norm2(psi1) ) << std::endl;
405 
406  D_w(psi1, chi, MINUS);
407  D_temp_prec2.unprecLinOp(psi2, chi, MINUS);
408  diff = psi2 - psi1;
409  QDPIO::cout << " D_2 - = " << sqrt( norm2(diff) / norm2(psi1) ) << std::endl;
410 
411  // ==============================================================================================
412  // ILU Prec Clover Op
413  // ===========================================================================================
414 
415  CloverFermActParams p;
416  p.Mass=0.01;
417  p.clovCoeffR = 0.91;
418  p.clovCoeffT = 1.07;
419  p.u0 = 1.0;
420 
421  p.anisoParam.anisoP = true;
422  p.anisoParam.t_dir = 3;
423  p.anisoParam.xi_0 = 2.464;
424  p.anisoParam.nu = 0.95;
425 
426 
427  QDPIO::cout << "ILU Prec Clover Op: " << std::endl;
428  ILUPrecSCprecTCloverLinOp D_temp_prec_clover(fs, p );
429 
430  UnprecCloverLinOp D_clov(fs, p);
431 
432  tmp1 = Gamma(15)*chi;
433  D_temp_prec_clover.invCRightLinOp(tmp2, tmp1, PLUS);
434  psi1 = Gamma(15)*tmp2;
435 
436  // psi2 = ( C_L^{-1} )^\dagger \chi
437  D_temp_prec_clover.invCLeftLinOp(psi2, chi, MINUS);
438  diff = psi2 - psi1;
439  QDPIO::cout << " Gamma5_1 = " << sqrt( norm2(diff) / norm2(psi1) ) << std::endl;
440 
441 
442  // Other way
443  // psi1 = \gamma_5 C_L^{-1} \gamma_5 \chi
444  tmp1 = Gamma(15)*chi;
445  D_temp_prec_clover.invCLeftLinOp(tmp2, tmp1, PLUS);
446  psi1 = Gamma(15)*tmp2;
447 
448  // psi2 = ( C_R^{-1} )^\dagger \chi
449  D_temp_prec_clover.invCRightLinOp(psi2, chi, MINUS);
450  diff = psi2 - psi1;
451  QDPIO::cout << " Gamma5_2 = " << sqrt( norm2(diff) / norm2(psi1) ) << std::endl;
452 
453  // Test LeftInverse is inverse of Left (Both orderings, PLUS and MINUS)
454  gaussian(chi);
455  D_temp_prec_clover.getFermBC().modifyF(chi);
456  D_temp_prec_clover.cLeftLinOp(psi1, chi, PLUS);
457  D_temp_prec_clover.invCLeftLinOp(tmp1, psi1, PLUS);
458 
459  diff = tmp1 - chi;
460  QDPIO::cout << " CL CL_INV + = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
461 
462 
463  D_temp_prec_clover.cLeftLinOp(psi1, chi, MINUS);
464  D_temp_prec_clover.invCLeftLinOp(tmp1, psi1, MINUS);
465 
466  diff = tmp1 - chi;
467  QDPIO::cout << " CL CL_INV - = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
468 
469  D_temp_prec_clover.invCLeftLinOp(psi1, chi, PLUS);
470  D_temp_prec_clover.cLeftLinOp(tmp1, psi1, PLUS);
471 
472  diff = tmp1 - chi;
473  QDPIO::cout << " CL_INV CL + = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
474 
475  D_temp_prec_clover.invCLeftLinOp(psi1, chi, MINUS);
476  D_temp_prec_clover.cLeftLinOp(tmp1, psi1, MINUS);
477 
478  diff = tmp1 - chi;
479  QDPIO::cout << " CL_INV CL - = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
480 
481 
482  // Test C_R is inverse of C_R_Inverse both orders, PLUS & Minus
483  gaussian(chi);
484  D_temp_prec_clover.getFermBC().modifyF(chi);
485  D_temp_prec_clover.cRightLinOp(psi1, chi, PLUS);
486  D_temp_prec_clover.invCRightLinOp(tmp1, psi1, PLUS);
487 
488  diff = tmp1 - chi;
489  QDPIO::cout << " CR CR_INV + = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
490 
491  D_temp_prec_clover.cRightLinOp(psi1, chi, MINUS);
492  D_temp_prec_clover.invCRightLinOp(tmp1, psi1, MINUS);
493 
494  diff = tmp1 - chi;
495  QDPIO::cout << " CR CR_INV - = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
496 
497  D_temp_prec_clover.invCRightLinOp(psi1, chi, PLUS);
498  D_temp_prec_clover.cRightLinOp(tmp1, psi1, PLUS);
499 
500  diff = tmp1 - chi;
501  QDPIO::cout << " CR_INV CR + = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
502 
503  D_temp_prec_clover.invCRightLinOp(psi1, chi, MINUS);
504  D_temp_prec_clover.cRightLinOp(tmp1, psi1, MINUS);
505 
506  diff = tmp1 - chi;
507  QDPIO::cout << " CR_INV CR - = " << sqrt(norm2(diff)/norm2(chi)) << std::endl;
508 
509 
510  // Now test against the normal D_op = C_L^{-1} unprecOp() C_R^{-1}
511  D_clov(psi1, chi, PLUS);
512  D_temp_prec_clover.unprecLinOp(psi2, chi, PLUS);
513  diff = psi2 - psi1;
514  QDPIO::cout << " D_2 + = " << sqrt( norm2(diff) / norm2(psi1) ) << std::endl;
515 
516  D_clov(psi1, chi, MINUS);
517  D_temp_prec_clover.unprecLinOp(psi2, chi, MINUS);
518  diff = psi2 - psi1;
519  QDPIO::cout << " D_2 - = " << sqrt( norm2(diff) / norm2(psi1) ) << std::endl;
520 
521 #if 0
522  // Schur Style preconditioning
523  EO3DPrecSCprecTWilsonLinOp D_schur_tprec(fs, Mass, aniso);
524 
525  QDPIO::cout << "Schur Style preconditioning tests " << std::endl;
526  QDPIO::cout << "================================= " << std::endl;
527 
528 
529 
530  for(int cb=0; cb < 2; cb++) {
531  gaussian(chi,rb3[cb]);
532 
533  tmp1[rb3[cb]] = Gamma(15)*chi;
534  D_schur_tprec.invCRightLinOp(tmp2, tmp1, PLUS,cb);
535  psi1[rb3[cb]] = Gamma(15)*tmp2;
536 
537  // psi2 = ( C_L^{-1} )^\dagger \chi
538  D_schur_tprec.invCLeftLinOp(psi2, chi, MINUS,cb);
539  diff[rb3[cb]] = psi2 - psi1;
540  QDPIO::cout << "cb="<<cb<<" Gamma5_1 = " << sqrt( norm2(diff,rb3[cb]) / norm2(psi1,rb3[cb]) ) << std::endl;
541 
542 
543  // Other way
544  // psi1 = \gamma_5 C_L^{-1} \gamma_5 \chi
545  tmp1[rb3[cb]] = Gamma(15)*chi;
546  D_schur_tprec.invCLeftLinOp(tmp2, tmp1, PLUS, cb);
547  psi1[rb3[cb]] = Gamma(15)*tmp2;
548 
549  // psi2 = ( C_R^{-1} )^\dagger \chi
550  D_schur_tprec.invCRightLinOp(psi2, chi, MINUS, cb);
551  diff[rb3[cb]] = psi2 - psi1;
552  QDPIO::cout << "cb="<<cb<< " Gamma5_2 = " << sqrt( norm2(diff, rb3[cb]) / norm2(psi1,rb3[cb]) ) << std::endl;
553 
554  // Test LeftInverse is inverse of Left (Both orderings, PLUS and MINUS)
555  gaussian(chi,rb3[cb]);
556  D_schur_tprec.cLeftLinOp(psi1, chi, PLUS,cb);
557  D_schur_tprec.invCLeftLinOp(tmp1, psi1, PLUS,cb);
558 
559  diff[rb3[cb]] = tmp1 - chi;
560  QDPIO::cout << "cb="<<cb<<" CL CL_INV + = " << sqrt(norm2(diff,rb3[cb])/norm2(chi,rb3[cb])) << std::endl;
561 
562 
563  D_schur_tprec.cLeftLinOp(psi1, chi, MINUS,cb);
564  D_schur_tprec.invCLeftLinOp(tmp1, psi1, MINUS,cb);
565 
566  diff[rb3[cb]] = tmp1 - chi;
567  QDPIO::cout << "cb="<<cb<<" CL CL_INV - = " << sqrt(norm2(diff,rb3[cb])/norm2(chi,rb3[cb])) << std::endl;
568 
569  D_schur_tprec.invCLeftLinOp(psi1, chi, PLUS,cb);
570  D_schur_tprec.cLeftLinOp(tmp1, psi1, PLUS,cb);
571 
572  diff[rb3[cb]] = tmp1 - chi;
573  QDPIO::cout << "cb="<<cb<<" CL_INV CL + = " << sqrt(norm2(diff,rb3[cb])/norm2(chi,rb3[cb])) << std::endl;
574 
575  D_schur_tprec.invCLeftLinOp(psi1, chi, MINUS,cb);
576  D_schur_tprec.cLeftLinOp(tmp1, psi1, MINUS,cb);
577 
578  diff[rb3[cb]] = tmp1 - chi;
579  QDPIO::cout << "cb="<<cb<<" CL_INV CL - = " << sqrt(norm2(diff,rb3[cb])/norm2(chi,rb3[cb])) << std::endl;
580 
581  // Test C_R is inverse of C_R_Inverse both orders, PLUS & Minus
582  gaussian(chi,rb3[cb]);
583  D_schur_tprec.cRightLinOp(psi1, chi, PLUS,cb);
584  D_schur_tprec.invCRightLinOp(tmp1, psi1, PLUS,cb);
585 
586  diff[rb3[cb]] = tmp1 - chi;
587  QDPIO::cout << "cb="<<cb<< " CR CR_INV + = " << sqrt(norm2(diff,rb3[cb])/norm2(chi,rb3[cb])) << std::endl;
588 
589  D_schur_tprec.cRightLinOp(psi1, chi, MINUS,cb);
590  D_schur_tprec.invCRightLinOp(tmp1, psi1, MINUS,cb);
591 
592  diff[rb3[cb]] = tmp1 - chi;
593  QDPIO::cout << "cb="<<cb<< " CR CR_INV - = " << sqrt(norm2(diff,rb3[cb])/norm2(chi,rb3[cb])) << std::endl;
594 
595  D_schur_tprec.invCRightLinOp(psi1, chi, PLUS,cb);
596  D_schur_tprec.cRightLinOp(tmp1, psi1, PLUS,cb);
597 
598  diff[rb3[cb]] = tmp1 - chi;
599  QDPIO::cout << "cb="<<cb<<" CR_INV CR + = " << sqrt(norm2(diff,rb3[cb])/norm2(chi,rb3[cb])) << std::endl;
600 
601  D_schur_tprec.invCRightLinOp(psi1, chi, MINUS,cb);
602  D_schur_tprec.cRightLinOp(tmp1, psi1, MINUS,cb);
603 
604  diff[rb3[cb]] = tmp1 - chi;
605  QDPIO::cout << "cb="<<cb<<" CR_INV CR - = " << sqrt(norm2(diff,rb3[cb])/norm2(chi,rb3[cb])) << std::endl;
606  }
607 
608  // D_w is unpreconditioned
609  // D_w1, chi, PLUS);
610  gaussian(chi);
611  D_schur_tprec.cRightLinOp(tmp1, chi, PLUS, 0);
612  D_schur_tprec.cRightLinOp(tmp1, chi, PLUS, 1);
613 
614  D_w(tmp2, tmp1, PLUS);
615 
616  D_schur_tprec.cLeftLinOp(psi1, tmp2, PLUS, 0);
617  D_schur_tprec.cLeftLinOp(psi1, tmp2, PLUS, 1);
618 
619  // --------------
620 
621  D_schur_tprec.evenEvenLinOp(psi2, chi, PLUS);
622  D_schur_tprec.evenOddLinOp(tmp1, chi, PLUS);
623  psi2[rb3[0]] += tmp1;
624 
625  D_schur_tprec.oddOddLinOp(psi2, chi, PLUS);
626  D_schur_tprec.oddEvenLinOp(tmp1, chi, PLUS);
627  psi2[rb3[1]] += tmp1;
628 
629 
630  for(int cb=0; cb < 2; cb++) {
631  diff[rb3[cb]] = psi2 - psi1;
632  QDPIO::cout << "cb="<<cb<<" D+ = " << sqrt(norm2(diff,rb3[cb])/norm2(chi,rb3[cb])) << std::endl;
633  }
634 
635  gaussian(chi);
636  D_w(psi1, chi, PLUS);
637  D_schur_tprec.unprecLinOp(psi2,chi,PLUS);
638  for(int cb=0; cb < 2; cb++) {
639  diff[rb3[cb]] = psi2 - psi1;
640  QDPIO::cout << "cb="<<cb<<" D+ = " << sqrt(norm2(diff,rb3[cb])/norm2(chi,rb3[cb])) << std::endl;
641  }
642 
643 
644  gaussian(chi);
645  D_schur_tprec.cLeftLinOp(tmp1, chi, MINUS, 0);
646  D_schur_tprec.cLeftLinOp(tmp1, chi, MINUS, 1);
647 
648  D_w(tmp2, tmp1, MINUS);
649 
650  D_schur_tprec.cRightLinOp(psi1, tmp2, MINUS, 0);
651  D_schur_tprec.cRightLinOp(psi1, tmp2, MINUS, 1);
652 
653  // --------------
654 
655 
656  D_schur_tprec.evenEvenLinOp(psi2, chi, MINUS);
657  D_schur_tprec.evenOddLinOp(tmp1, chi, MINUS);
658  psi2[rb3[0]] += tmp1;
659 
660  D_schur_tprec.oddOddLinOp(psi2, chi, MINUS);
661  D_schur_tprec.oddEvenLinOp(tmp1, chi, MINUS);
662  psi2[rb3[1]] += tmp1;
663 
664 
665  for(int cb=0; cb < 2; cb++) {
666  diff[rb3[cb]] = psi2 - psi1;
667  QDPIO::cout << "cb="<<cb<<" D- = " << sqrt(norm2(diff,rb3[cb])/norm2(chi,rb3[cb])) << std::endl;
668  }
669 
670  gaussian(chi);
671  D_w(psi1, chi, MINUS);
672  D_schur_tprec.unprecLinOp(psi2,chi,MINUS);
673  for(int cb=0; cb < 2; cb++) {
674  diff[rb3[cb]] = psi2 - psi1;
675  QDPIO::cout << "cb="<<cb<<" D- = " << sqrt(norm2(diff,rb3[cb])/norm2(chi,rb3[cb])) << std::endl;
676  }
677 
678 
679 
680  // -----------------------------------------------------------------
681  EO3DPrecSCprecTCloverLinOp D_schur_clov_tprec(fs,p);
682  gaussian(chi);
683  D_clov(psi1, chi, PLUS);
684  D_schur_clov_tprec.unprecLinOp(psi2,chi,PLUS);
685  for(int cb=0; cb < 2; cb++) {
686  diff[rb3[cb]] = psi2 - psi1;
687  QDPIO::cout << "cb="<<cb<<" D+ = " << sqrt(norm2(diff,rb3[cb])/norm2(chi,rb3[cb])) << std::endl;
688  }
689 
690  gaussian(chi);
691  D_clov(psi1, chi, MINUS);
692  D_schur_clov_tprec.unprecLinOp(psi2,chi,MINUS);
693  for(int cb=0; cb < 2; cb++) {
694  diff[rb3[cb]] = psi2 - psi1;
695  QDPIO::cout << "cb="<<cb<<" D- = " << sqrt(norm2(diff,rb3[cb])/norm2(chi,rb3[cb])) << std::endl;
696  }
697 #endif
698 
699  // Time to bolt
700  Chroma::finalize();
701 
702 #endif
703 
704  exit(0);
705 }
chroma.h
Primary include file for CHROMA in application codes.
Chroma::FermionAction::createState
virtual FermState< T, P, Q > * createState(const Q &q) const
Given links (coordinates Q) create the state needed for the linear operators.
Definition: fermact.h:59
Chroma::Handle
Class for counted reference semantics.
Definition: handle.h:33
Chroma::ILUPrecSpaceCentralPrecTimeWilsonTypeFermAct
ILUPreconditioned Spatial, Central Temporal Preconditioned Wilson-like fermion actions with derivativ...
Definition: iluprec_s_cprec_t_wilstype_fermact_w.h:19
Chroma::ILUPrecSpaceCentralPrecTimeWilsonTypeFermAct::linOp
virtual ILUPrecSpaceCentralPrecTimeLinearOperator< T, P, Q > * linOp(Handle< FermState< T, P, Q > > state) const =0
Produce a linear operator for this action.
Chroma::SingletonHolder::Instance
static T & Instance()
Definition: singleton.h:432
Chroma::UnprecCloverLinOp
Unpreconditioned Clover-Dirac operator.
Definition: unprec_clover_linop_w.h:25
Chroma::UnprecWilsonFermAct
Unpreconditioned Wilson fermion action.
Definition: unprec_wilson_fermact_w.h:32
Chroma::UnprecWilsonFermAct::linOp
UnprecLinearOperator< T, P, Q > * linOp(Handle< FermState< T, P, Q > > state) const
Produce a linear operator for this action.
Definition: unprec_wilson_fermact_w.cc:70
Chroma::UnprecWilsonLinOp
Unpreconditioned Wilson-Dirac operator.
Definition: unprec_wilson_linop_w.h:41
eo3dprec_s_cprec_t_clover_linop_w.h
eo3dprec_s_cprec_t_wilson_linop_w.h
Chroma::read
void read(XMLReader &xml, const std::string &path, AsqtadFermActParams &param)
Read parameters.
Definition: asqtad_fermact_params_s.cc:33
Chroma::unitarityCheck
void unitarityCheck(const multi1d< LatticeColorMatrixF3 > &u)
Check the unitarity of color matrix in SU(N)
Definition: unit_check.cc:20
Chroma::HotSt
void HotSt(multi1d< LatticeColorMatrix > &u)
Set a gauge field from a sample of (almost) Haar measure.
Definition: hotst.cc:34
Chroma::gaugeStartup
void gaugeStartup(XMLReader &gauge_file_xml, XMLReader &gauge_xml, multi1d< LatticeColorMatrix > &u, Cfg_t &cfg)
Initialize the gauge fields.
Definition: gauge_startup.cc:36
Chroma::proginfo
void proginfo(XMLWriter &xml)
Print out basic information about this program.
Definition: proginfo.cc:24
Chroma::CFG_TYPE_WEAK_FIELD
@ CFG_TYPE_WEAK_FIELD
Definition: enum_cfgtype_io.h:36
Nd
Nd
Definition: meslate.cc:74
tmp2
Double tmp2
Definition: mesq.cc:30
Chroma::InlineDiscoEigCGEnv::P
multi1d< LatticeColorMatrix > P
Definition: inline_disco_eigcg_w.cc:315
Chroma::WilsonTypeFermActsEnv::registerAll
bool registerAll()
Register all the factories.
Definition: fermacts_aggregate_w.cc:196
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::Mass
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > & Mass
Definition: pbg5p_w.cc:29
Chroma::push
push(xml_out,"Condensates")
Chroma::T
LinOpSysSolverMGProtoClover::T T
Definition: syssolver_linop_clover_mg_proto.cc:63
Chroma::i
int i
Definition: pbg5p_w.cc:55
Chroma::initialize
void initialize(int *argc, char ***argv)
Chroma initialisation routine.
Definition: chroma_init.cc:114
Chroma::MINUS
@ MINUS
Definition: chromabase.h:45
Chroma::PLUS
@ PLUS
Definition: chromabase.h:45
Chroma::finalize
void finalize(void)
Chroma finalization routine.
Definition: chroma_init.cc:308
Chroma::chi
multi1d< LatticeFermion > chi(Ncb)
Chroma::pop
pop(xml_out)
Chroma::MesPlq
void MesPlq(const multi1d< LatticeColorMatrixF3 > &u, multi2d< Double > &plane_plaq, Double &link)
Definition: mesplq.cc:83
Chroma::cb
int cb
Definition: invbicg.cc:120
Chroma::getXMLOutputInstance
XMLFileWriter & getXMLOutputInstance()
Get xml output instance.
Definition: chroma_init.cc:359
Chroma::zero
Double zero
Definition: invbicg.cc:106
testing::internal::string
::std::string string
Definition: gtest.h:1979
Chroma::Cfg_t
Gauge configuration structure.
Definition: cfgtype_io.h:16
Chroma::CloverFermActParams
Params for clover ferm acts.
Definition: clover_fermact_params_w.h:16
linkage_hack
bool linkage_hack()
To insure linking of code, place the registered code flags here.
Definition: t_temp_prec.cc:17
main
int main(int argc, char **argv)
Definition: t_temp_prec.cc:27
Generated by doxygen 1.9.1