S_w->refresh(gauge_state);
S_w->refreshInternalFields(gauge_state);
int main(int argc, char *argv[]) {  // Initialise QDP  QDP_initialize(&argc, &argv);  // Setup a small lattice  const int nrow_arr[] = {2, 2, 2, 2};  multi1d<int> nrow(Nd);  nrow=nrow_arr;  Layout::setLattSize(nrow);  Layout::create();  multi1d<LatticeColorMatrix> u(Nd);  {    XMLReader file_xml;    XMLReader config_xml;    Cfg_t foo; foo.cfg_type=CFG_TYPE_DISORDERED;    // Cfg_t foo; foo.cfg_type=CFG_TYPE_SZIN; foo.cfg_file="./CFGIN";    gaugeStartup(file_xml, config_xml, u, foo);  }  // Dump output  XMLFileWriter xml_out("./XMLDAT");  push(xml_out, "t_monomial");  // Read Parameters  std::string monomial_name;           // String for Factory  XMLReader param_in("DATA");  Handle< ExactMonomial<multi1d<LatticeColorMatrix>,     multi1d<LatticeColorMatrix> > > S_w;  try {         // Snarf it all    XMLReader paramtop(param_in, "/MonomialTest");    // Get the string for the factory    read(paramtop, "Monomial", S_w);  }  catch(const string& e) {     QDPIO::cerr << "Error reading XML" << endl;    QDP_abort(1);  }  // Fictitious momenta for now  multi1d<LatticeColorMatrix> p(Nd);  for(int mu=0; mu<Nd; mu++) {     gaussian(p[mu]);  }  // Create a field state  GaugeFieldState gauge_state(p,u);  // Refresh Pseudofermions  S_w->refresh(gauge_state);  // Compute Force from Monomial  multi1d<LatticeColorMatrix> dsdq(Nd);  S_w->dsdq(dsdq, gauge_state);  // Compute action from monomial  Double monomial_S = S_w->S(gauge_state);  QDPIO::cout << "monomial_S = " << monomial_S << endl;  // Check against normal version  // Downcast for direct handling  UnprecTwoFlavorWilsonFermMonomial& S_down = dynamic_cast<UnprecTwoFlavorWilsonFermMonomial&>(*S_w);  // Use the debug accessors  const LatticeFermion& phi = S_down.debugGetPhi();  LatticeFermion X=zero;  // This will repeat the solve... need to put caching in at some stage  S_down.debugGetX(X, gauge_state);  // My energy measurement  Double my_S = innerProductReal(phi,X);  QDPIO::cout << "My S = " << my_S << endl;  // Compute force the old fashioned way  // Get at the FermAct  const UnprecWilsonFermAct& S_f = S_down.debugGetFermAct();  multi1d<LatticeColorMatrix> dsdq2(Nd);    // Call the old force routine  wilson_dsdu(S_f, dsdq2, gauge_state.getQ(), X);  // Compare to new force   for(int mu=0; mu < Nd; mu++) {     taproj(dsdq[mu]);    taproj(dsdq2[mu]);        push(xml_out, "dsdu");    write(xml_out, "dsdu_1", dsdq[mu]);    write(xml_out, "dsdu_2", dsdq2[mu]);    pop(xml_out);        LatticeColorMatrix dsdu_diff=dsdq[mu] - dsdq2[mu];        Double sum_diff=norm2(dsdu_diff);    QDPIO::cout << "Mu = " << mu << " Sum Diff=" << sum_diff << endl;        push(xml_out, "ForceDiff");    write(xml_out, "mu", mu);    write(xml_out, "dsdu_diff", dsdu_diff);    pop(xml_out);  }    // End   pop(xml_out);  xml_out.close();    // Finish  QDP_finalize();  exit(0);}