CHROMA
t_minvert_qphix.cc
Go to the documentation of this file.
1 // $Id: t_invert4_precwilson.cc,v 3.4 2009-10-09 13:59:46 bjoo Exp $
2 
3 #include "chroma.h"
4 #include <iostream>
5 #include <sstream>
6 #include <iomanip>
7 #include <string>
8 
9 #include <cstdio>
10 
11 #include <stdlib.h>
12 #include <sys/time.h>
13 #include <math.h>
14 
15 
16 #include "io/xml_group_reader.h"
19 #include <string>
20 
24 //#include "actions/ferm/invert/multi_syssolver_mdagm_cg_chrono_clover.h"
26 using namespace Chroma;
27 
28 struct AppParams {
29  multi1d<int> nrow;
30  Cfg_t inputCfg;
31  GroupXML_t fermact;
32  // MultiSysSolverCGChronoCloverParams invParam;
34 };
35 
36 
37 
38 void checkInverter(const AppParams& p, multi1d<LatticeColorMatrix>& u)
39 {
40 
41  typedef LatticeFermion T;
42  typedef multi1d<LatticeColorMatrix> Q;
43  typedef multi1d<LatticeColorMatrix> P;
44 
45  typedef LatticeFermionF TF;
46  typedef multi1d<LatticeColorMatrixF> QF;
47  typedef multi1d<LatticeColorMatrixF> PF;
48 
49  std::istringstream is(p.fermact.xml);
50  XMLReader fermact_xml(is);
51  Handle< WilsonTypeFermAct<T,P,Q> > S_f ( TheWilsonTypeFermActFactory::Instance().createObject(p.fermact.id, fermact_xml, p.fermact.path) );
52 
53 
54  //S_f(
55  // new EvenOddPrecCloverFermAct(CreateFermStateEnv::reader(fermact_xml, p.fermact.path), CloverFermActParams(fermact_xml, p.fermact.path)));
56 
57 
58  Handle< FermState<T,P,Q> > connect_state( S_f->createState(u) );
59  Handle< LinearOperator<LatticeFermion> > D_op( S_f->linOp(connect_state) );
60 
61 
62  // MdagMMultiSysSolverCGChronoClover mprec(D_op, connect_state, p.invParam);
63  MdagMMultiSysSolverQPhiXClover<T,LatticeColorMatrix> mprec(D_op, connect_state, p.invParam);
64  int n_shifts=12;
65  multi1d<Real> shifts(n_shifts);
66  shifts[0] = 1.51480322230288e-05 ;
67  shifts[1] = 0.000165259114300968 ;
68  shifts[2] = 0.000661184012887753 ;
69  shifts[3] = 0.00214955953261573 ;
70  shifts[4] = 0.00657190811374048 ;
71  shifts[5] = 0.0197043312297316 ;
72  shifts[6] = 0.0587640674448055 ;
73  shifts[7] = 0.175534749714768 ;
74  shifts[8] = 0.529999019144126 ;
75  shifts[9] = 1.6565878795396 ;
76  shifts[10] = 5.78532483906374 ;
77  shifts[11] = 30.6835527589452 ;
78 
79  multi1d<Real> RsdCG(shifts.size());
80  if (p.invParam.RsdTarget.size() == 1) {
81  RsdCG = p.invParam.RsdTarget[0];
82  }
83  else if (p.invParam.RsdTarget.size() == RsdCG.size()) {
84 
85  RsdCG = p.invParam.RsdTarget;
86  }
87  else {
88 
89  QDPIO::cerr << "MdagMMultiSysSolverCGChronoClover: shifts incompatible" << std::endl;
90  QDP_abort(1);
91  }
92 
93 
94 
95  LatticeFermion chi;
96  gaussian(chi, D_op->subset());
97 
98  // Initialized psi
99  multi1d<LatticeFermion> psi(n_shifts);
100  for(int i=0; i < n_shifts; i++){
101  psi[i] = zero;
102  }
103 
105  StopWatch swatch;
106  swatch.reset();
107  swatch.start();
108  MInvCG2(*D_op, chi, psi, shifts, RsdCG, p.invParam.MaxIter, res.n_count);
109  swatch.stop();
110 
111  // Check solutions
112  for(int i=0; i < n_shifts; i++) {
113  LatticeFermion r;
114  r[D_op->subset()]= chi;
115  LatticeFermion tmp1,tmp2;
116  (*D_op)(tmp1, psi[i], PLUS);
117  (*D_op)(tmp2, tmp1, MINUS);
118  tmp2[D_op->subset()] += shifts[i]*psi[i];
119  r[D_op->subset()] -= tmp2;
120  Double norm = norm2(r, D_op->subset())/norm2(chi, D_op->subset());
121 
122  QDPIO::cout << "soln "<< i <<" : RsdCG=" << RsdCG[i] << " r=" << sqrt(norm) << std::endl;
123  }
124  QDPIO::cout << "MinvCG: n_count = " << res.n_count << " Time=" << swatch.getTimeInSeconds() << std::endl;
125 
126  swatch.reset();
127  swatch.start();
128 
129  SystemSolverResults_t res2 = mprec(psi, shifts, chi);
130 
131  swatch.stop();
132  for(int i=0; i < n_shifts; i++) {
133  LatticeFermion r;
134  r[D_op->subset()]= chi;
135  LatticeFermion tmp1,tmp2;
136  (*D_op)(tmp1, psi[i], PLUS);
137  (*D_op)(tmp2, tmp1, MINUS);
138  tmp2[D_op->subset()] += shifts[i]*psi[i];
139  r[D_op->subset()] -= tmp2;
140  Double norm = norm2(r, D_op->subset())/norm2(chi, D_op->subset());
141 
142  QDPIO::cout << "soln "<< i <<" : RsdCG=" << RsdCG[i] << " r=" << sqrt(norm) << std::endl;
143  }
144  QDPIO::cout << "Chrono: n_count = " << res2.n_count << " Time=" << swatch.getTimeInSeconds() << std::endl;
145 
146 }
147 
148 
149 void read(XMLReader& r, const std::string path, AppParams& p)
150 {
151  read(r, "nrow", p.nrow);
152  read(r, "Cfg", p.inputCfg);
153  p.fermact = readXMLGroup(r, "FermionAction", "FermAct");
154  read(r, "InvertParam", p.invParam);
155 
156 }
157 
158 bool linkageHack(void)
159 {
160  bool foo = true;
161 
162  // Inline Measurements
164  foo &= GaugeInitEnv::registerAll();
165 
166  return foo;
167 }
168 
169 int main(int argc, char **argv)
170 {
171  // Put the machine into a known state
172  Chroma::initialize(&argc, &argv);
173  QDPIO::cout << "Linkage = " << linkageHack() << std::endl;
174 
175 
177 
178  XMLReader xml_in;
179  try
180  {
181  xml_in.open(Chroma::getXMLInputFileName());
182  read(xml_in, "/Param", params);
183  }
184  catch(const std::string& e)
185  {
186  QDPIO::cerr << "Caught Exception reading XML: " << e << std::endl;
187  QDP_abort(1);
188  }
189  Layout::setLattSize(params.nrow);
190  Layout::create();
191 
192  multi1d<LatticeColorMatrix> u(Nd);
193  XMLReader gauge_file_xml, gauge_xml;
194  gaugeStartup(gauge_file_xml, gauge_xml, u, params.inputCfg);
195  unitarityCheck(u);
196 
197  // Setup the lattice
198 
199  XMLFileWriter& xml_out = Chroma::getXMLOutputInstance();
200  push(xml_out,"t_invert");
201  push(xml_out,"params");
202  write(xml_out, "nrow", params.nrow);
203  pop(xml_out); // Params
204 
205  // Measure the plaquette on the gauge
206  MesPlq(xml_out, "Observables", u);
207  xml_out.flush();
208 
210 
211  pop(xml_out);
212  xml_out.close();
213 
215 
216  exit(0);
217 }
Primary include file for CHROMA in application codes.
Class for counted reference semantics.
Definition: handle.h:33
static T & Instance()
Definition: singleton.h:432
Parameters for Clover fermion action.
Even-odd preconditioned Clover fermion linear operator.
All ferm create-state method.
void read(XMLReader &xml, const std::string &path, AsqtadFermActParams &param)
Read parameters.
void write(XMLWriter &xml, const std::string &path, const AsqtadFermActParams &param)
Writer parameters.
void unitarityCheck(const multi1d< LatticeColorMatrixF3 > &u)
Check the unitarity of color matrix in SU(N)
Definition: unit_check.cc:20
void gaugeStartup(XMLReader &gauge_file_xml, XMLReader &gauge_xml, multi1d< LatticeColorMatrix > &u, Cfg_t &cfg)
Initialize the gauge fields.
void MInvCG2(const LinearOperator< LatticeFermionF > &M, const LatticeFermionF &chi, multi1d< LatticeFermionF > &psi, const multi1d< RealF > &shifts, const multi1d< RealF > &RsdCG, int MaxCG, int &n_count)
Definition: minvcg2.cc:376
GroupXML_t readXMLGroup(XMLReader &xml_in, const std::string &path, const std::string &type_name)
Read group and return as a std::string.
Params params
Nd
Definition: meslate.cc:74
Double tmp2
Definition: mesq.cc:30
Multishift Conjugate-Gradient algorithm for a Linear Operator.
Solve a M*psi=chi linear system by BiCGStab.
bool registerAll()
Register all the factories.
bool registerAll()
Register all the factories.
multi1d< LatticeColorMatrix > P
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
gaussian(aux)
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > enum InvType invType const multi1d< Real > & RsdCG
Definition: pbg5p_w.cc:30
static multi1d< LatticeColorMatrix > u
LinOpSysSolverMGProtoClover::Q Q
push(xml_out,"Condensates")
LinOpSysSolverMGProtoClover::T T
int i
Definition: pbg5p_w.cc:55
void initialize(int *argc, char ***argv)
Chroma initialisation routine.
Definition: chroma_init.cc:114
@ MINUS
Definition: chromabase.h:45
@ PLUS
Definition: chromabase.h:45
void finalize(void)
Chroma finalization routine.
Definition: chroma_init.cc:308
multi1d< LatticeFermion > chi(Ncb)
bool linkageHack(void)
Definition: const_hmc.cc:660
std::string getXMLInputFileName()
Get input file name.
Definition: chroma_init.cc:88
LatticeFermion psi
Definition: mespbg5p_w.cc:35
pop(xml_out)
void MesPlq(const multi1d< LatticeColorMatrixF3 > &u, multi2d< Double > &plane_plaq, Double &link)
Definition: mesplq.cc:83
XMLFileWriter & getXMLOutputInstance()
Get xml output instance.
Definition: chroma_init.cc:359
const WilsonTypeFermAct< multi1d< LatticeFermion > > & S_f
Definition: pbg5p_w.cc:27
Double zero
Definition: invbicg.cc:106
FloatingPoint< double > Double
Definition: gtest.h:7351
::std::string string
Definition: gtest.h:1979
int norm
Definition: qtopcor.cc:35
MultiSysSolverQPhiXCloverParams invParam
Gauge configuration structure.
Definition: cfgtype_io.h:16
Hold group xml and type id.
Holds return info from SystemSolver call.
Definition: syssolver.h:17
Factory for solving M*psi=chi where M is not hermitian or pos. def.
std::string fermact_xml
void checkInverter(const AppParams &p, multi1d< LatticeColorMatrix > &u)
int main(int argc, char **argv)
multi1d< LatticeColorMatrixF > QF
Definition: t_quda_tprec.cc:19
multi1d< LatticeColorMatrixF > PF
Definition: t_quda_tprec.cc:20
LatticeFermionF TF
Definition: t_quda_tprec.cc:17
Read an XML group as a std::string.