multi2d<DComplex> pion(16, input.param,nrow[3]);
multi2d<DComplex> pion(16, input.param.nrow[3]);
int main(int argc, char **argv){  // Put the machine into a known state  QDP_initialize(&argc, &argv);  // Input parameter structure  Propagator_input_t  input;  // Instantiate xml reader for DATA  XMLReader xml_in("DATA");  // Read data  read(xml_in, "/propagator", input);  // Specify lattice size, shape, etc.  Layout::setLattSize(input.param.nrow);  Layout::create();  // Read in the configuration along with relevant information.  multi1d<LatticeColorMatrix> u(Nd);    XMLReader gauge_xml;  switch (input.param.cfg_type)   {  case CFG_TYPE_NERSC :    readArchiv(gauge_xml, u, input.cfg.cfg_file);    break;  default :    QDP_error_exit("Configuration type is unsupported.");  }  //   //  gauge invariance test  //    // gauge transformed gauge fields//  multi1d<LatticeColorMatrix> u_trans(Nd);  // gauge transform//  LatticeColorMatrix v ;  //  gaussian(v);//  reunit(v) ; //  for(int dir = 0 ; dir < Nd ; ++dir)//    {//      u_trans[dir] = v*u[dir]*adj(shift(v,FORWARD,dir)) ;//    }  // Read in the source along with relevant information.  LatticePropagator quark_prop_source;  XMLReader source_xml;  switch (input.param.prop_type)   {  case PROP_TYPE_SZIN ://    readSzinQprop(source_xml, quark_prop_source, input.prop.source_file);    quark_prop_source = zero;    break;  default :    QDP_error_exit("Propagator type is unsupported.");  }  // Instantiate XML writer for XMLDAT  XMLFileWriter xml_out("XMLDAT");  push(xml_out, "propagator");  // Write out the input  write(xml_out, "Input", xml_in);  // Write out the config header  write(xml_out, "Config_info", gauge_xml);  // Write out the source header  write(xml_out, "Source_info", source_xml);  push(xml_out, "Output_version");  write(xml_out, "out_version", 1);  pop(xml_out);  xml_out.flush();  // Check if the gauge field configuration is unitarized  unitarityCheck(u);  // Calculate some gauge invariant observables just for info.  Double w_plaq, s_plaq, t_plaq, link;  MesPlq(u, w_plaq, s_plaq, t_plaq, link);  push(xml_out, "Observables");  Write(xml_out, w_plaq);  Write(xml_out, s_plaq);  Write(xml_out, t_plaq);  Write(xml_out, link);  pop(xml_out);  // Fix to the coulomb gauge  int n_gf;  int j_decay = Nd-1;  coulGauge(u, n_gf, j_decay, input.param.GFAccu, input.param.GFMax, true, input.param.OrPara);  QDPIO::cout << "No. of gauge fixing iterations =" << n_gf << endl;  // Calcluate plaq on the gauge fixed field  MesPlq(u, w_plaq, s_plaq, t_plaq, link);  push(xml_out, "Is this gauge invariant?");  Write(xml_out, w_plaq);  Write(xml_out, s_plaq);  Write(xml_out, t_plaq);  Write(xml_out, link);  pop(xml_out);  xml_out.flush();  // Create a fermion BC. Note, the handle is on an ABSTRACT type.  Handle< FermBC<LatticeFermion> >  fbc(new SimpleFermBC<LatticeFermion>(input.param.boundary));  //  // Initialize fermion action  //  AsqtadFermAct S_f(fbc, input.param.Mass, input.param.u0);  // Set up a state for the current u,  // (compute fat & triple links)  // Use S_f.createState so that S_f can pass in u0  Handle<const ConnectState > state(S_f.createState(u));  Handle<const EvenOddLinearOperator<LatticeFermion> > D_asqtad(S_f.linOp(state));  // Create a fermion to apply linop to.  LatticeFermion tmp1, tmp2;  tmp1 = zero;//  Test walfil code//  for(int src_ind = 0; src_ind < 8; ++src_ind){//    walfil(tmp1, 0, 3, 0, src_ind);//      tmp2  =  zero;      // Apply Linop//      (*D_asqtad).evenOddLinOp(tmp2, tmp1, PLUS); ////       push(xml_out, "dslash");//       Write(xml_out, tmp1);//       Write(xml_out, tmp2);//       pop(xml_out);//    }   Handle<const LinearOperator<LatticeFermion> > MdagM_asqtad(S_f.lMdagM(state));  //  // Loop over the source color, creating the source  // and calling the relevant propagator routines. The QDP  // terminology is that a staggered propagator is a matrix in color space  //   //  LatticePropagator quark_propagator;  XMLBufferWriter xml_buf;  int ncg_had = 0;  int n_count;  LatticeFermion q_source, psi;  multi1d<LatticePropagator> stag_prop(8);//  for(int t_source = 0; t_source < 17; t_source += 16) {    for(int src_ind = 0; src_ind < 8; ++src_ind){      psi = zero;   // note this is ``zero'' and not 0      int t_source = 0;      for(int color_source = 0; color_source < Nc; ++color_source) {        QDPIO::cout << "Inversion for Color =  " << color_source << endl;        QDPIO::cout << "Source time slice = " << t_source << endl;        int spin_source = 0;        q_source = zero ;//        srcfil(q_source, t_src_even, color_source, 0);//  Use a wall source        walfil(q_source, t_source, j_decay, color_source, src_ind);        // Use the last initial guess as the current guess        // Compute the propagator for given source color/spin         // int n_count;        S_f.qprop(psi, state, q_source, CG_INVERTER,                   input.param.RsdCG, input.param.MaxCG, n_count);            ncg_had += n_count;              push(xml_out,"Qprop");        write(xml_out, "Mass" , input.param.Mass);        write(xml_out, "RsdCG", input.param.RsdCG);        Write(xml_out, n_count);        pop(xml_out);        /*         * Move the solution to the appropriate components         * of quark propagator.        */        FermToProp(psi, quark_propagator, color_source, spin_source);      }  //color_source          stag_prop[src_ind] = quark_propagator;      } // end src_ind        multi2d<DComplex> pion(16, input.param,nrow[3]);      staggeredPionsFollana(stag_prop, pion, j_decay);    push(xml_out, "Here are all 16 pions");      for(int i=0; i < NUM_STAG_PIONS; i++) {      ostringstream tag;      tag << "pion" << i;      push(xml_out, tag.str());      Write(xml_out, pion[i]);      pop(xml_out);      }    // Instantiate XML buffer to make the propagator header    XMLBufferWriter prop_xml;    push(prop_xml, "propagator");    // Write out the input    write(prop_xml, "Input", xml_in);    // Write out the config header    write(prop_xml, "Config_info", gauge_xml);    // Write out the source header    write(prop_xml, "Source_info", source_xml);    pop(prop_xml);  // Save the propagator//    switch (input.param.prop_type) //    {//    case PROP_TYPE_SZIN://     writeSzinQprop(quark_propagator, input.prop.prop_file, input.param.Mass);//    break;//  case PROP_TYPE_SCIDAC://    writeQprop(prop_xml, quark_propagator, input.prop.prop_file);//    break; //    default ://     QDP_error_exit("Propagator type is unsupported.");//    }//  } //t_source;  xml_out.close();  xml_in.close();  // Time to bolt  QDP_finalize();  exit(0);}