Layout::setLattSize(params.nrow);
Layout::setLattSize(input.param.nrow);
int main(int argc, char **argv){  // Put the machine into a known state  QDP_initialize(&argc, &argv);  int G5 = Ns*Ns - 1;  // Read the parameters   Param_t params;  try {     readParams("./DATA", params);  }  catch(const string& s) {     QDPIO::cerr << "Caught exception " << s << endl;    exit(1);  }  // Setup the lattice  Layout::setLattSize(params.nrow);  Layout::create();  // Write out the params  XMLFileWriter xml_out("t_ovlap5d.xml");  push(xml_out, "overlapTest");  dumpParams(xml_out, params);    // The Gauge Field  multi1d<LatticeColorMatrix> u(Nd);    switch ((GaugeStartType)params.gauge_start_type) {   case COLD_START:    for(int j = 0; j < Nd; j++) {       u(j) = Real(1);    }    break;  case HOT_START:    // Hot (disordered) start    for(int j=0; j < Nd; j++) {       random(u(j));      reunit(u(j));    }    break;  case FILE_START:    // Start from File     switch( (GaugeFormat)params.gauge_file_format) {     case SZIN_GAUGE_FORMAT:      {	XMLReader szin_xml;	readSzin(szin_xml, u, params.gauge_filename);	try { 	  push(xml_out, "GaugeInfo");	  xml_out << szin_xml;	  pop(xml_out);	}	catch(const string& e) {	  cerr << "Error: " << e << endl;	}	      }      break;    case NERSC_GAUGE_FORMAT:      {	XMLReader nersc_xml;	readArchiv(nersc_xml, u, params.gauge_filename);	try { 	  push(xml_out, "GaugeInfo");	  xml_out << nersc_xml;	  pop(xml_out);	}	catch(const string& e) {	  cerr << "Error: " << e << endl;	}      }      break;    default:      ostringstream file_read_error;      file_read_error << "Unknown gauge file format" << params.gauge_file_format ;      throw file_read_error.str();    }    break;  default:    ostringstream startup_error;    startup_error << "Unknown start type " << params.gauge_start_type <<endl;    throw startup_error.str();  }  // Measure the plaquette on the gauge  Double w_plaq, s_plaq, t_plaq, link;  MesPlq(u, w_plaq, s_plaq, t_plaq, link);  push(xml_out, "plaquette");  write(xml_out, "w_plaq", w_plaq);  write(xml_out, "s_plaq", s_plaq);  write(xml_out, "t_plaq", t_plaq);  write(xml_out, "link", link);  pop(xml_out);  //! Wilsoniums;  // Put this puppy into a handle to allow Zolo to copy it around as a **BASE** class  // WARNING: the handle now owns the data. The use of a bare S_w below is legal,  // but just don't delete it.  // Create a FermBC  Handle< FermBC<LatticeFermion> >  fbc(new SimpleFermBC<LatticeFermion>(params.boundary));   Handle< FermBC< multi1d<LatticeFermion> > >  fbc5(new SimpleFermBC<multi1d<LatticeFermion> >(params.boundary));    Handle<UnprecWilsonTypeFermAct<LatticeFermion> >  S_w(new UnprecWilsonFermAct(fbc, params.wilson_mass));    XMLBufferWriter my_writer;  //! N order Zolo approx, with wilson action.  Zolotarev5DFermActArray   S(fbc5,			      S_w, 			      params.quark_mass,			      params.approx_order, 			      my_writer);     const ConnectState* connect_state_ptr;  multi1d<LatticeFermion> eigen_vecs;  if( params.lambda.size() == 0 ) {     // Connect State with no eigenvectors    connect_state_ptr = S.createState(u, 				      Real(params.approx_min), 				      Real(params.approx_max));  }  else {    // Connect State with eigenvectors    if( params.szin_eig ) {      readEigenVecs(u, 		    dynamic_cast<UnprecWilsonFermAct&>(*S_w), 		    params.lambda, 		    eigen_vecs, 		    params.wilson_mass, 		    params.szin_eig, 		    xml_out);    }    else {      QDP_error_exit("Non SZIN e-values not yet implmeneted");    }    connect_state_ptr = S.createState(u, 				      params.lambda,				      eigen_vecs,				      params.lambda_max);  }    // Stuff the pointer into a handle. Now, the handle owns the data.  Handle<const ConnectState> connect_state(connect_state_ptr);    // Make me a linop (this callls the initialise function)  Handle<const LinearOperator< multi1d< LatticeFermion > > > D_op(S.linOp(connect_state));  push(xml_out, "Zolotarev5DInternal");  xml_out << my_writer;  pop(xml_out);  int N5 = S.size();  LatticeFermion chi4;  LatticeFermion psi4;  LatticeFermion tmp;  int n_count;  gaussian(chi4);  chi4 /= sqrt(norm2(chi4));   // We solve M_5d psi = gamma_5 chi  multi1d<LatticeFermion> psi( N5 );  multi1d<LatticeFermion> chi( N5 );  for(int i = 0; i < N5; i++) {     psi[i] = zero;    chi[i] = zero;  }   chi[N5-1] = Gamma(G5)*chi4;             // 4D source in last component                                         // is gamma5 chi  multi1d<LatticeFermion> tmp5_1( N5 );  // CGNE  //  //  tmp5_1 = ( M^{dag} M )^{-1} chi  //  //  M^{dag} M is hermitian so M^{dag}M = M M^{dag}  //  // hence tmp5_1 = M^{-dag} M^{-1} chi    // Put solution into psi  check inverse  (*D_op)(tmp5_1, chi, MINUS);  InvCG2( *D_op, tmp5_1, psi, params.rsd_cg, params.max_cg, n_count);   // Multiply back to check inverse  (*D_op)(tmp5_1, psi, PLUS);    multi1d<LatticeFermion> tmp5_2( N5 );  Double dnorm = Double(0);  Double g5chi_norm = Double(0);  for(int i = 0; i < N5; i++) {     tmp5_2[i] = chi[i] - tmp5_1[i];    dnorm += norm2(tmp5_2[i]);    g5chi_norm += norm2(chi[i]);  }      Double r_norm5 = sqrt(dnorm/g5chi_norm);  QDPIO::cout << "|| chi - D psi ||/ || chi || = " << r_norm5 << endl;  push(xml_out, "Inv5DCheck");  write(xml_out, "r_norm", r_norm5);  pop(xml_out);  // Get the bottom component into psi4  psi4 = psi[N5-1];  // Scale to get right normalsation.   // Note this only works because we gamma_5-d the source initially  psi4 *= Real(2)/(Real(1)-params.quark_mass);    // This is the residuum in 4D  LatticeFermion r;  //! N order Zolo approx, with wilson action.  Zolotarev4DFermAct   S4(fbc,			  S_w, 			  params.quark_mass,			  params.approx_order, 			  1.0e-9,			  500,			  my_writer);  Handle<const LinearOperator< LatticeFermion > > D_op4(S4.linOp(connect_state));    (*D_op4)(r, psi4, PLUS);  r -= chi4;             // Ungamma5-ed  Double r_norm4 = sqrt(norm2(r)/norm2(chi4));  QDPIO::cout << "|| chi4 - D psi4 || = "<< r_norm4 << endl;  push(xml_out, "Inv4Dcheck");  write(xml_out, "r_norm", r_norm4);  pop(xml_out); #if 0   // Qprop Test  multi1d<int> coord(Nd);  coord[0]=0; coord[1] = 0; coord[2] = 0; coord[3] = 0;  chi4 = zero;  tmp = zero;  // Point source  srcfil(chi4, coord, 0, 0);  psi4 = zero;    S4.qprop(psi4, connect_state, chi4, CG_INVERTER, params.rsd_cg, params.max_cg, n_count);  // Check original solution. psi4 contains 1/(1-m)[D^{-1} - 1 ]  psi4 *= (Real(1)-params.quark_mass);  psi4 += chi4; // add back contact term  psi4 = D^{-1} chi4  (*D_op4)(r, psi4, PLUS);         // r = D psi4  r -= chi4;                       // r = D psi4 - chi4  r_norm4 = sqrt(norm2(r)/norm2(chi4));  QDPIO::cout << "4D Qprop || r || = " << r_norm4 << endl;  push(xml_out, "Consistency4D");  write(xml_out, "r_norm", r_norm4);  pop(xml_out);  LatticeFermion psi4_2=zero;  S.qprop(psi4_2, connect_state, chi4, CG_INVERTER, params.rsd_cg, params.max_cg, n_count);  psi4_2 *= (Real(1)-params.quark_mass);  // Multiply out overall normalisation  psi4_2 += chi4;                         // add back in contact term  (*D_op4)(r, psi4_2, PLUS);         // r = D psi4  r -= chi4;                       // r = D psi4 - chi4  r_norm4 = sqrt(norm2(r)/norm2(chi4));   QDPIO::cout << "5D Qprop || r || = " << r_norm4 << endl;  push(xml_out, "Consistency5D");  write(xml_out, "r_norm", r_norm4);  pop(xml_out);  // Remove contact_term  psi4 -= chi4;  psi4_2 -= chi4;  //  Overlall normalisation  psi4 *= Real(1)/( Real(1) - params.quark_mass);  psi4_2 *= Real(1)/( Real(1) - params.quark_mass);  r = psi4 - psi4_2;  r_norm4 = sqrt(norm2(r)/norm2(chi4));  QDPIO::cout << "|| 4D Qprop - 5D Qprop || = " << r_norm4 << endl;  push(xml_out, "Qprop4D5Ddiff");  write(xml_out, "r_norm", r_norm4);  pop(xml_out);#endif  pop(xml_out);  QDP_finalize();      exit(0);}