#include <config.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include <util/qcdio.h>
#ifdef PARALLEL
#include <comms/sysfunc_cps.h>
#endif
#include <comms/scu.h>
#include <comms/glb.h>

#include <util/lattice.h>
#include <util/time_cps.h>
#include <util/smalloc.h>

#include <util/command_line.h>

#include<unistd.h>
#include<config.h>

#include <alg/lanc_arg.h>
#include <util/spincolorflavormatrix.h>
#include <alg/propmanager.h>
#include <alg/fix_gauge_arg.h>
#include <alg/alg_gparitycontract.h>
#include <alg/prop_dft.h>

#include <map>
#include <string>
#include <sstream>
#include <bitset>
#include <algorithm>

#include <util/omp_wrapper.h>

CPS_START_NAMESPACE

void PropDFT::global_coord(const int &site, int *into_vec){    
  int rem = site;
  for(int i=0;i<4;i++){
    into_vec[i] = rem % GJP.NodeSites(i) + GJP.NodeCoor(i)*GJP.NodeSites(i);
    rem /= GJP.NodeSites(i);
  }
}

//Add a sink momentum to the set of those generated by the Fourier transform
//Automatically adds minus the vector as this allows for optimization in retrieval when the bilinear desired is the complex or Hermitian conjugate of one already calculated
void PropDFT::add_momentum(std::vector<Float> sink_mom){
  mom_idx_map[ sink_mom ] = nmom++;
  if(sink_mom[0]!=0.0 || sink_mom[1]!=0.0 || sink_mom[2]!=0.0){
    sink_mom[0]*=-1; sink_mom[1]*=-1; sink_mom[2]*=-1;
    mom_idx_map[ sink_mom ] = nmom++;
  }
}

int PropDFT::momIdxMinusP(const int &pidx) const{
  std::vector<Float> mom;
  bool found(false);
  for(mom_idx_map_type::const_iterator mom_it = mom_idx_map.begin(); mom_it != mom_idx_map.end(); mom_it++){
    if(mom_it->second == pidx){
      mom = mom_it->first;
      found = true;
    }
  }
  if(!found) ERR.General("PropDFT","momIdxMinusP","Could not find momentum with idx %d\n",pidx);

  for(int i=0;i<3;i++) mom[i]*=-1;
    
  mom_idx_map_type::const_iterator conj_it = mom_idx_map.find(mom);
  if(conj_it == mom_idx_map.end()) ERR.General("PropDFT","momIdxMinusP","Minus mom not present in map, why?");
  return conj_it->second;
}



static bool mom_sort_pred(const std::pair<Float,int> &i, const std::pair<Float,int> &j){ 
  return (i.first<j.first);
}

void PropDFT::find_p2sorted(std::vector< std::pair<Float,int> > &p2list, std::map<int,std::vector<Float> > &p2map){
  //Find all p^2 and sort, output into p2list and p2map
  p2list.clear(); p2map.clear(); p2list.reserve(nmom);
  for(mom_idx_map_type::iterator mom_it = mom_idx_map.begin(); mom_it != mom_idx_map.end(); ++mom_it){
    const std::vector<Float> & mom = mom_it->first;
    const int &vec_pos = mom_it->second;
    std::pair<Float,int> p2idx( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2], vec_pos );
    p2list.push_back(p2idx);
    p2map[vec_pos] = mom;
  }
  std::sort(p2list.begin(),p2list.end(),mom_sort_pred);
}



void _FourierProp_helper<SpinColorFlavorMatrix>::site_matrix(SpinColorFlavorMatrix &into, QPropWcontainer &prop, Lattice &lat, const int &site){
  into.generate(prop, lat,site);
}

void _FourierProp_helper<SpinColorFlavorMatrix>::lattice_sum(SpinColorFlavorMatrix &what){
  _FourierProp_helper<WilsonMatrix>::lattice_sum(what(0,0));
  _FourierProp_helper<WilsonMatrix>::lattice_sum(what(0,1));
  _FourierProp_helper<WilsonMatrix>::lattice_sum(what(1,0));
  _FourierProp_helper<WilsonMatrix>::lattice_sum(what(1,1));
}

void _FourierProp_helper<SpinColorFlavorMatrix>::mult_gauge_fix_mat(SpinColorFlavorMatrix &what, const int &site, Lattice &lat){
  Matrix const* gfmat_f0 = lat.FixGaugeMatrix(site,0);
  Matrix const* gfmat_f1 = lat.FixGaugeMatrix(site,1);
  if(gfmat_f0 == NULL || gfmat_f1 == NULL) ERR.General("_FourierProp_helper<SpinColorFlavorMatrix>","mult_gauge_fix_mat(SpinColorFlavorMatrix &what, const int &site, Lattice &lat)","No gauge fixing matrix for site %d",site);
  what(0,0).LeftTimesEqual(*gfmat_f0);
  what(0,1).LeftTimesEqual(*gfmat_f0);
  what(1,0).LeftTimesEqual(*gfmat_f1);
  what(1,1).LeftTimesEqual(*gfmat_f1);
}

void _FourierProp_helper<SpinColorFlavorMatrix>::write(FILE *fp, const SpinColorFlavorMatrix &mat, const Float &p2, const std::vector<Float> &mom, const int &t){
  //format is p^2 px py pz t s0 c0 f0 s1 c1 f1 val.re val.im
  for(int s0=0;s0<4;s0++) for(int c0=0;c0<3;c0++) for(int f0=0;f0<2;f0++)
  for(int s1=0;s1<4;s1++) for(int c1=0;c1<3;c1++) for(int f1=0;f1<2;f1++){
    const Complex &val = mat(s0,c0,f0,s1,c1,f1);
    Fprintf(fp, "%.16e %.16e %.16e %.16e %d %d %d %d %d %d %d %.16e %.16e\n",p2,mom[0],mom[1],mom[2],t,s0,c0,f0,s1,c1,f1,val.real(),val.imag());
  }
}

void _FourierProp_helper<WilsonMatrix>::site_matrix(WilsonMatrix &into, QPropWcontainer &prop, Lattice &lat, const int &site){
  into = prop.getProp(lat).SiteMatrix(site,0);
}

void _FourierProp_helper<WilsonMatrix>::lattice_sum(WilsonMatrix &what){
  Float* w = (Float*)what.ptr(); //returns Rcomplex*
  static const int size = 2*12*12;
  slice_sum(w, size, 99); //99 is a *magic* number (we are abusing slice_sum here)
}

void _FourierProp_helper<WilsonMatrix>::mult_gauge_fix_mat(WilsonMatrix &what, const int &site, Lattice &lat){
  Matrix const* gfmat = lat.FixGaugeMatrix(site,0);
  if(gfmat == NULL) ERR.General("_FourierProp_helper<WilsonMatrix>","mult_gauge_fix_mat(WilsonMatrix &what, const int &site, Lattice &lat)","No gauge fixing matrix for site %d",site);
  what.LeftTimesEqual(*gfmat);
}

void _FourierProp_helper<WilsonMatrix>::write(FILE *fp, const WilsonMatrix &mat, const Float &p2, const std::vector<Float> &mom, const int &t){
  //format is p^2 px py pz t s0 c0 s1 c1 val.re val.im
  for(int s0=0;s0<4;s0++) for(int c0=0;c0<3;c0++)
  for(int s1=0;s1<4;s1++) for(int c1=0;c1<3;c1++){
    const Complex &val = mat(s0,c0,s1,c1);
    Fprintf(fp, "%.16e %.16e %.16e %.16e %d %d %d %d %d %.16e %.16e\n",p2,mom[0],mom[1],mom[2],t,s0,c0,s1,c1,val.real(),val.imag());
  }
}

const int _PropagatorBilinear_helper<SpinColorFlavorMatrix>::nidx(64); //16*4

//currently only allow for spin and flavour matrices between the propagators
int _PropagatorBilinear_helper<SpinColorFlavorMatrix>::scf_map(const int &spinidx, const int &flavidx){ return flavidx + 4*spinidx; }

//unmap a spin-flavour index into separate spin and flavour indices
std::pair<int,int> _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unmap(const int &scf_idx){
  std::pair<int,int> out; int rem = scf_idx;
  out.second = rem % 4; rem/=4; //flav
  out.first = rem; //spin
  return out;
}
//Calculate the coefficient of the result of transposing and/or conjugating the spin-flavor matrix
Float _PropagatorBilinear_helper<SpinColorFlavorMatrix>::coeff(const int &scf_idx, const bool &transpose, const bool &conj){ 
  std::pair<int,int> sfidx = unmap(scf_idx);
  return AlgGparityContract::qdp_gcoeff(sfidx.first,transpose,conj) * AlgGparityContract::pauli_coeff(sfidx.second,transpose,conj);
}

void _PropagatorBilinear_helper<SpinColorFlavorMatrix>::site_matrix(SpinColorFlavorMatrix &into, QPropWcontainer &prop, Lattice &lat, const int &site){ 
  return _FourierProp_helper<SpinColorFlavorMatrix>::site_matrix(into,prop,lat,site); 
}

//right-multiply with spin and flavour matrices
void _PropagatorBilinear_helper<SpinColorFlavorMatrix>::rmult_matrix(SpinColorFlavorMatrix &into, const std::pair<int,int> &spin_flav){
  const static FlavorMatrixType fmap[4] = {sigma0, sigma1, sigma2, sigma3};   
  AlgGparityContract::qdp_gr(into,spin_flav.first);      
  into.pr(fmap[spin_flav.second]);
}

void _PropagatorBilinear_helper<SpinColorFlavorMatrix>::lattice_sum(SpinColorFlavorMatrix &what){
  return _FourierProp_helper<SpinColorFlavorMatrix>::lattice_sum(what);
}
  
void _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unit_matrix(SpinColorFlavorMatrix &mat){
  mat = 0.0;
  for(int f=0;f<2;f++) for(int s=0;s<4;s++) for(int c=0;c<3;c++) mat(s,c,f,s,c,f) = 1.0;
}

void _PropagatorBilinear_helper<SpinColorFlavorMatrix>::write(FILE *fp, const SpinColorFlavorMatrix &mat, const int &idx, const Float &p2, const std::vector<Float> &mom, const int &t){
  std::pair<int,int> gamma_sigma = unmap(idx);
  
  //format is gamma sigma p^2 px py pz t s0 c0 f0 s1 c1 f1 val.re val.im
  for(int s0=0;s0<4;s0++) for(int c0=0;c0<3;c0++) for(int f0=0;f0<2;f0++)
  for(int s1=0;s1<4;s1++) for(int c1=0;c1<3;c1++) for(int f1=0;f1<2;f1++){
    const Complex &val = mat(s0,c0,f0,s1,c1,f1);
    Fprintf(fp, "%d %d %.16e %.16e %.16e %.16e %d %d %d %d %d %d %d %.16e %.16e\n",gamma_sigma.first,gamma_sigma.second,p2,mom[0],mom[1],mom[2],t,s0,c0,f0,s1,c1,f1,val.real(),val.imag());
  }	
}

const int _PropagatorBilinear_helper<WilsonMatrix>::nidx(16);

int _PropagatorBilinear_helper<WilsonMatrix>::scf_map(const int &spinidx, const int &flavidx){ return spinidx; }

//put spin index on first element of output
std::pair<int,int> _PropagatorBilinear_helper<WilsonMatrix>::unmap(const int &sc_idx){
  return std::pair<int,int>(sc_idx,0);
}

//Calculate the coefficient of the result of transposing and/or conjugating the spin matrix
Float _PropagatorBilinear_helper<WilsonMatrix>::coeff(const int &sidx, const bool &transpose, const bool &conj){ 
  return AlgGparityContract::qdp_gcoeff(sidx,transpose,conj); 
}
typedef _PropagatorBilinear_generics::map_info_scmat map_info_type; //the map type

void _PropagatorBilinear_helper<WilsonMatrix>::site_matrix(WilsonMatrix &into, QPropWcontainer &prop, Lattice &lat, const int &site){ 
  return _FourierProp_helper<WilsonMatrix>::site_matrix(into,prop,lat,site); 
}
  
//right-multiply with spin matrices
void _PropagatorBilinear_helper<WilsonMatrix>::rmult_matrix(WilsonMatrix &into, const std::pair<int,int> &spin_flav){   
  AlgGparityContract::qdp_gr(into,spin_flav.first);      
}

void _PropagatorBilinear_helper<WilsonMatrix>::lattice_sum(WilsonMatrix &what){
  return _FourierProp_helper<WilsonMatrix>::lattice_sum(what);
}
void _PropagatorBilinear_helper<WilsonMatrix>::unit_matrix(WilsonMatrix &mat){
  mat = 0.0;
  for(int s=0;s<4;s++) for(int c=0;c<3;c++) mat(s,c,s,c) = 1.0;
}

void _PropagatorBilinear_helper<WilsonMatrix>::write(FILE *fp, const WilsonMatrix &mat, const int &idx, const Float &p2, const std::vector<Float> &mom, const int &t){
  //format is gamma p^2 px py pz t s0 c0 s1 c1 val.re val.im
  for(int s0=0;s0<4;s0++) for(int c0=0;c0<3;c0++)
  for(int s1=0;s1<4;s1++) for(int c1=0;c1<3;c1++){
    const Complex &val = mat(s0,c0,s1,c1);
    Fprintf(fp, "%d %.16e %.16e %.16e %.16e %d %d %d %d %d %.16e %.16e\n",idx,p2,mom[0],mom[1],mom[2],t,s0,c0,s1,c1,val.real(),val.imag());
  }	
}


void _ContractedBilinear_helper<SpinColorFlavorMatrix>::write(FILE *fp, const Rcomplex &val, const int &idx1, const int &idx2, const Float &p2, const std::vector<Float> &mom, const int &t){
  std::pair<int,int> gamma_sigma1 = _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unmap(idx1);
  std::pair<int,int> gamma_sigma2 = _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unmap(idx2);

  //format is gamma1 sigma1 gamma2 sigma2 p^2 px py pz t val.re val.im
  Fprintf(fp, "%d %d %d %d %.16e %.16e %.16e %.16e %d %.16e %.16e\n",
	  gamma_sigma1.first,gamma_sigma1.second,
	  gamma_sigma2.first,gamma_sigma2.second,
	  p2,mom[0],mom[1],mom[2],t,val.real(),val.imag());
}


void _ContractedBilinear_helper<WilsonMatrix>::write(FILE *fp, const Rcomplex &val, const int &idx1, const int &idx2, const Float &p2, const std::vector<Float> &mom, const int &t){
  //format is gamma1 gamma2 p^2 px py pz t val.re val.im
  Fprintf(fp, "%d %d %.16e %.16e %.16e %.16e %d %.16e %.16e\n",idx1,idx2,p2,mom[0],mom[1],mom[2],t,val.real(),val.imag());
}


void _ContractedBilinear_helper<SpinColorFlavorMatrix>::binary_write(FILE *fp, const Rcomplex &val, const int &idx1, const int &idx2, const Float &p2, const std::vector<Float> &mom, const int &t){
  std::pair<int,int> gamma_sigma1 = _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unmap(idx1);
  std::pair<int,int> gamma_sigma2 = _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unmap(idx2);

  if(!UniqueID()){
    int ibuf[5] = { gamma_sigma1.first, gamma_sigma1.second, gamma_sigma2.first, gamma_sigma2.second, t };
    fwrite(ibuf, sizeof(int), 5, fp);
    double dbuf[6] = { p2,mom[0],mom[1],mom[2],val.real(),val.imag() };
    fwrite(dbuf, sizeof(double), 6 ,fp);
  }
}

void _ContractedBilinear_helper<WilsonMatrix>::binary_write(FILE *fp, const Rcomplex &val, const int &idx1, const int &idx2, const Float &p2, const std::vector<Float> &mom, const int &t){
  if(!UniqueID()){
    int ibuf[3] = { idx1,idx2, t };
    fwrite(ibuf, sizeof(int), 3, fp);
    double dbuf[6] = { p2,mom[0],mom[1],mom[2],val.real(),val.imag() };
    fwrite(dbuf, sizeof(double), 6 ,fp);
  }
}





void _ContractedBilinear_helper<SpinColorFlavorMatrix>::write(std::ostream &fp, const Rcomplex &val, const int &idx1, const int &idx2, const Float &p2, const std::vector<Float> &mom, const int &t){
  std::pair<int,int> gamma_sigma1 = _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unmap(idx1);
  std::pair<int,int> gamma_sigma2 = _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unmap(idx2);

  //format is gamma1 sigma1 gamma2 sigma2 p^2 px py pz t val.re val.im
  if(!UniqueID())
    fp << gamma_sigma1.first<< " " << gamma_sigma1.second<< " " << gamma_sigma2.first<< " " << gamma_sigma2.second<< " " 
       << p2<< " " << mom[0]<< " " << mom[1]<< " " << mom[2]<< " " << t<< " " << val.real()<< " " << val.imag() << std::endl;
}


void _ContractedBilinear_helper<WilsonMatrix>::write(std::ostream &fp, const Rcomplex &val, const int &idx1, const int &idx2, const Float &p2, const std::vector<Float> &mom, const int &t){
  //format is gamma1 gamma2 p^2 px py pz t val.re val.im
  if(!UniqueID())
    fp << idx1<< " " << idx2<< " " << p2<< " " << mom[0]<< " " << mom[1]<< " " << mom[2]<< " " << t<< " " << val.real()<< " " << val.imag() << std::endl;
}




//24^4 tensor in spin,color and flavor
int QuadrilinearSCFVertex::map(const int &s1,const int &c1, const int &f1,
			       const int &s2,const int &c2, const int &f2,
			       const int &s3,const int &c3, const int &f3,
			       const int &s4,const int &c4, const int &f4) const{
  return s1 + 4*(c1 + 3*(f1 + 2*(
				 s2 + 4*(c2 + 3*(f2 + 2*(
							 s3 + 4*(c3 + 3*(f3 + 2*(				   
										 s4 + 4*(c4 + 3*(f4
												 )))))))))));
}
void QuadrilinearSCFVertex::unmap(int idx,
				  int &s1,int &c1, int &f1,
				  int &s2,int &c2, int &f2,
				  int &s3,int &c3, int &f3,
				  int &s4,int &c4, int &f4) const{
  s1 = idx%4; idx/=4;
  c1 = idx%3; idx/=3;
  f1 = idx%2; idx/=2;
  s2 = idx%4; idx/=4;
  c2 = idx%3; idx/=3;
  f2 = idx%2; idx/=2;
  s3 = idx%4; idx/=4;
  c3 = idx%3; idx/=3;
  f3 = idx%2; idx/=2;
  s4 = idx%4; idx/=4;
  c4 = idx%3; idx/=3;
  f4 = idx;
}



QuadrilinearSCFVertex::QuadrilinearSCFVertex(): size(24*24*24*24), tensor(new Rcomplex[24*24*24*24]){}

QuadrilinearSCFVertex::QuadrilinearSCFVertex(const Rcomplex &val): size(24*24*24*24), tensor(new Rcomplex[24*24*24*24]){
  for(int i=0;i<size;i++) tensor[i] = val;
}

  //form tensor from outer product of 2 SpinColorFlavorMatrix
QuadrilinearSCFVertex::QuadrilinearSCFVertex(const SpinColorFlavorMatrix &bilA, const SpinColorFlavorMatrix &bilB, const Rcomplex &phase): 
  size(24*24*24*24), tensor(new Rcomplex[24*24*24*24]){
  bool phase_mult(false); if(phase.real()!=1.0 || phase.imag()!=0.0) phase_mult = true;

  for(int s1=0; s1<4; ++s1) for(int c1=0; c1<3; ++c1) for(int f1=0; f1<2; ++f1)
  for(int s2=0; s2<4; ++s2) for(int c2=0; c2<3; ++c2) for(int f2=0; f2<2; ++f2)
  for(int s3=0; s3<4; ++s3) for(int c3=0; c3<3; ++c3) for(int f3=0; f3<2; ++f3)
  for(int s4=0; s4<4; ++s4) for(int c4=0; c4<3; ++c4) for(int f4=0; f4<2; ++f4){
    int idx = map(s1,c1,f1,s2,c2,f2,s3,c3,f3,s4,c4,f4);
    tensor[idx] = bilA(s1,c1,f1,s2,c2,f2)*bilB(s3,c3,f3,s4,c4,f4);
    if(phase_mult) tensor[idx] *= phase;
  }
}
  //add tensor formed from outer product of 2 SpinColorFlavorMatrix
QuadrilinearSCFVertex &QuadrilinearSCFVertex::add(const SpinColorFlavorMatrix &bilA, const SpinColorFlavorMatrix &bilB, const Rcomplex &phase){
  bool phase_mult(false); if(phase.real()!=1.0 || phase.imag()!=0.0) phase_mult = true;

  for(int s1=0; s1<4; ++s1) for(int c1=0; c1<3; ++c1) for(int f1=0; f1<2; ++f1)
  for(int s2=0; s2<4; ++s2) for(int c2=0; c2<3; ++c2) for(int f2=0; f2<2; ++f2)
  for(int s3=0; s3<4; ++s3) for(int c3=0; c3<3; ++c3) for(int f3=0; f3<2; ++f3)
  for(int s4=0; s4<4; ++s4) for(int c4=0; c4<3; ++c4) for(int f4=0; f4<2; ++f4){
    int idx = map(s1,c1,f1,s2,c2,f2,s3,c3,f3,s4,c4,f4);
    Rcomplex toadd = bilA(s1,c1,f1,s2,c2,f2)*bilB(s3,c3,f3,s4,c4,f4);
    if(phase_mult) toadd *= phase;
    tensor[idx] += toadd;
  }
  return *this;
}

QuadrilinearSCFVertex& QuadrilinearSCFVertex::operator=(const Rcomplex &val){
  for(int i=0;i<size;i++) tensor[i] = val;
  return *this;
}
QuadrilinearSCFVertex& QuadrilinearSCFVertex::operator=(const Float &val){
  for(int i=0;i<size;i++) tensor[i] = val;
  return *this;
}

QuadrilinearSCFVertex& QuadrilinearSCFVertex::operator+=(const QuadrilinearSCFVertex &rhs){
  for(int i=0;i<size;i++) tensor[i]+=rhs.tensor[i];
  return *this;
}

const Rcomplex & QuadrilinearSCFVertex::operator()(const int &s1,const int &c1, const int &f1,
						   const int &s2,const int &c2, const int &f2,
						   const int &s3,const int &c3, const int &f3,
						   const int &s4,const int &c4, const int &f4) const{
  return tensor[ map(s1,c1,f1,s2,c2,f2,s3,c3,f3,s4,c4,f4) ];
}

Rcomplex* QuadrilinearSCFVertex::ptr(){ return tensor; }

QuadrilinearSCFVertex::~QuadrilinearSCFVertex(){ delete[] tensor; }



//16^4 tensor in spin and color

int QuadrilinearSCVertex::map(const int &s1,const int &c1,
			      const int &s2,const int &c2,
			      const int &s3,const int &c3,
			      const int &s4,const int &c4) const{
  return s1 + 4*(c1 + 3*(
			 s2 + 4*(c2 + 3*(
					 s3 + 4*(c3 + 3*(
							 s4 + 4*(c4
								 )))))));
		   
  ;
}
void QuadrilinearSCVertex::unmap(int idx,
				 int &s1,int &c1,
				 int &s2,int &c2,
				 int &s3,int &c3,
				 int &s4,int &c4) const{
  s1 = idx%4; idx/=4;
  c1 = idx%3; idx/=3;
  s2 = idx%4; idx/=4;
  c2 = idx%3; idx/=3;
  s3 = idx%4; idx/=4;
  c3 = idx%3; idx/=3;
  s4 = idx%4; idx/=4;
  c4 = idx;
}


QuadrilinearSCVertex::QuadrilinearSCVertex(): size(12*12*12*12), tensor(new Rcomplex[12*12*12*12]){}
QuadrilinearSCVertex::QuadrilinearSCVertex(const Rcomplex &val): size(12*12*12*12), tensor(new Rcomplex[12*12*12*12]){
  for(int i=0;i<size;i++) tensor[i] = val;
}

  //form tensor from outer product of 2 WilsonMatrix
QuadrilinearSCVertex::QuadrilinearSCVertex(const WilsonMatrix &bilA, const WilsonMatrix &bilB, const Rcomplex &phase): 
  size(12*12*12*12), tensor(new Rcomplex[12*12*12*12]){
  bool phase_mult(false); if(phase.real()!=1.0 || phase.imag()!=0.0) phase_mult = true;

  for(int s1=0; s1<4; ++s1) for(int c1=0; c1<3; ++c1)
  for(int s2=0; s2<4; ++s2) for(int c2=0; c2<3; ++c2)
  for(int s3=0; s3<4; ++s3) for(int c3=0; c3<3; ++c3)
  for(int s4=0; s4<4; ++s4) for(int c4=0; c4<3; ++c4){
    int idx = map(s1,c1,s2,c2,s3,c3,s4,c4);
    tensor[idx] = bilA(s1,c1,s2,c2)*bilB(s3,c3,s4,c4);
    if(phase_mult) tensor[idx] *= phase;
  }
}

//add tensor formed from outer product of 2 WilsonMatrix
QuadrilinearSCVertex &QuadrilinearSCVertex::add(const WilsonMatrix &bilA, const WilsonMatrix &bilB, const Rcomplex &phase){
  bool phase_mult(false); if(phase.real()!=1.0 || phase.imag()!=0.0) phase_mult = true;

  for(int s1=0; s1<4; ++s1) for(int c1=0; c1<3; ++c1)
  for(int s2=0; s2<4; ++s2) for(int c2=0; c2<3; ++c2)
  for(int s3=0; s3<4; ++s3) for(int c3=0; c3<3; ++c3)
  for(int s4=0; s4<4; ++s4) for(int c4=0; c4<3; ++c4){
    int idx = map(s1,c1,s2,c2,s3,c3,s4,c4);
    Rcomplex toadd = bilA(s1,c1,s2,c2)*bilB(s3,c3,s4,c4);
    if(phase_mult) toadd *= phase;
    tensor[idx] += toadd;
  }
  return *this;
}

QuadrilinearSCVertex& QuadrilinearSCVertex::operator=(const Rcomplex &val){
  for(int i=0;i<size;i++) tensor[i] = val;
  return *this;
}
QuadrilinearSCVertex& QuadrilinearSCVertex::operator=(const Float &val){
  for(int i=0;i<size;i++) tensor[i] = val;
  return *this;
}

QuadrilinearSCVertex& QuadrilinearSCVertex::operator+=(const QuadrilinearSCVertex &rhs){
  for(int i=0;i<size;i++) tensor[i]+=rhs.tensor[i];
  return *this;
}

const Rcomplex & QuadrilinearSCVertex::operator()(const int &s1,const int &c1,
						  const int &s2,const int &c2,
						  const int &s3,const int &c3,
						  const int &s4,const int &c4) const{
  return tensor[ map(s1,c1,s2,c2,s3,c3,s4,c4) ];
}

Rcomplex* QuadrilinearSCVertex::ptr(){ return tensor; }

QuadrilinearSCVertex::~QuadrilinearSCVertex(){ delete[] tensor; }





void _PropagatorQuadrilinear_helper<SpinColorFlavorMatrix>::write(FILE *fp, const QuadrilinearSCFVertex &tensor,  const int &idx1, const int &idx2, const Float &p2, const std::vector<Float> &mom, const int &t){
  std::pair<int,int> gamma_sigma1 = _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unmap(idx1);
  std::pair<int,int> gamma_sigma2 = _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unmap(idx2);

  //format is gamma1 sigma1 gamma2 sigma2 p^2 px py pz t s0 c0 f0 s1 c1 f1 s2 c2 f2 s3 c3 f3 val.re val.im
  for(int s0=0;s0<4;s0++) for(int c0=0;c0<3;c0++) for(int f0=0;f0<2;f0++)
  for(int s1=0;s1<4;s1++) for(int c1=0;c1<3;c1++) for(int f1=0;f1<2;f1++)
  for(int s2=0;s2<4;s2++) for(int c2=0;c2<3;c2++) for(int f2=0;f2<2;f2++)
  for(int s3=0;s3<4;s3++) for(int c3=0;c3<3;c3++) for(int f3=0;f3<2;f3++){
    const Complex &val = tensor(s0,c0,f0,s1,c1,f1,s2,c2,f2,s3,c3,f3);
    Fprintf(fp, "%d %d %d %d %.16e %.16e %.16e %.16e %d %d %d %d %d %d %d %d %d %d %d %d %d %.16e %.16e\n",
	    gamma_sigma1.first,gamma_sigma1.second,
	    gamma_sigma2.first,gamma_sigma2.second,
	    p2,mom[0],mom[1],mom[2],t,s0,c0,f0,s1,c1,f1,s2,c2,f2,s3,c3,f3,val.real(),val.imag());
  }
}



void _PropagatorQuadrilinear_helper<SpinColorFlavorMatrix>::lattice_sum(QuadrilinearSCFVertex &t){
  Float* w = (Float*)t.ptr(); //returns Rcomplex*
  static const int size = 2*24*24*24*24;
  slice_sum(w, size, 99); //99 is a *magic* number (we are abusing slice_sum here)
}


void _PropagatorQuadrilinear_helper<WilsonMatrix>::write(FILE *fp, const QuadrilinearSCVertex &tensor, const int &idx1, const int &idx2, const Float &p2, const std::vector<Float> &mom, const int &t){
  //format is gamma1 gamma2 p^2 px py pz t s0 c0 s1 c1 s2 c2 s3 c3 val.re val.im
  for(int s0=0;s0<4;s0++) for(int c0=0;c0<3;c0++)
  for(int s1=0;s1<4;s1++) for(int c1=0;c1<3;c1++)
  for(int s2=0;s2<4;s2++) for(int c2=0;c2<3;c2++)
  for(int s3=0;s3<4;s3++) for(int c3=0;c3<3;c3++){
    const Complex &val = tensor(s0,c0,s1,c1,s2,c2,s3,c3);
    Fprintf(fp, "%d %d %.16e %.16e %.16e %.16e %d %d %d %d %d %d %d %d %d %.16e %.16e\n",
	    idx1,idx2,
	    p2,mom[0],mom[1],mom[2],t,s0,c0,s1,c1,s2,c2,s3,c3,val.real(),val.imag());
  }
}



void _PropagatorQuadrilinear_helper<WilsonMatrix>::lattice_sum(QuadrilinearSCVertex &t){
  Float* w = (Float*)t.ptr(); //returns Rcomplex*
  static const int size = 2*12*12*12*12;
  slice_sum(w, size, 99); //99 is a *magic* number (we are abusing slice_sum here)
}


void _ContractedWallSinkBilinearSpecMomentum_helper<SpinColorFlavorMatrix>::write(FILE *fp, const Rcomplex &val, const int &idx1, const int &idx2, const Float &p2, 
										  const std::vector<Float> &mom1, const std::vector<Float> &mom2, const int &t){
  std::pair<int,int> gamma_sigma1 = _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unmap(idx1);
  std::pair<int,int> gamma_sigma2 = _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unmap(idx2);

  //format is gamma1 sigma1 gamma2 sigma2 p^2 p1x p1y p1z p2x p2y p2z t val.re val.im
  Fprintf(fp, "%d %d %d %d %.16e %.16e %.16e %.16e %.16e %.16e %.16e %d %.16e %.16e\n",
	  gamma_sigma1.first,gamma_sigma1.second,
	  gamma_sigma2.first,gamma_sigma2.second,
	  p2,mom1[0],mom1[1],mom1[2],mom2[0],mom2[1],mom2[2],t,val.real(),val.imag());
}


void _ContractedWallSinkBilinearSpecMomentum_helper<WilsonMatrix>::write(FILE *fp, const Rcomplex &val, const int &idx1, const int &idx2, const Float &p2, 
									 const std::vector<Float> &mom1, const std::vector<Float> &mom2, const int &t){
  //format is gamma1 gamma2 p^2 p1x p1y p1z p2x p2y p2z t val.re val.im
  Fprintf(fp, "%d %d %.16e %.16e %.16e %.16e %.16e %.16e %.16e %d %.16e %.16e\n",idx1,idx2,p2,mom1[0],mom1[1],mom1[2],mom2[0],mom2[1],mom2[2],t,val.real(),val.imag());
}



void _ContractedWallSinkBilinearSpecMomentum_helper<SpinColorFlavorMatrix>::binary_write(FILE *fp, const Rcomplex &val, const int &idx1, const int &idx2, const Float &p2, 
										  const std::vector<Float> &mom1, const std::vector<Float> &mom2, const int &t){
  std::pair<int,int> gamma_sigma1 = _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unmap(idx1);
  std::pair<int,int> gamma_sigma2 = _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unmap(idx2);

  if(!UniqueID()){
    int ibuf[5] = { gamma_sigma1.first, gamma_sigma1.second, gamma_sigma2.first, gamma_sigma2.second, t };
    fwrite(ibuf, sizeof(int), 5, fp);
    double dbuf[9] = { p2,mom1[0],mom1[1],mom1[2],mom2[0],mom2[1],mom2[2],val.real(),val.imag() };
    fwrite(dbuf, sizeof(double), 9 ,fp);
  }
}


void _ContractedWallSinkBilinearSpecMomentum_helper<WilsonMatrix>::binary_write(FILE *fp, const Rcomplex &val, const int &idx1, const int &idx2, const Float &p2, 
									 const std::vector<Float> &mom1, const std::vector<Float> &mom2, const int &t){
  if(!UniqueID()){
    int ibuf[3] = { idx1,idx2,t };
    fwrite(ibuf, sizeof(int), 3, fp);
    double dbuf[9] = { p2,mom1[0],mom1[1],mom1[2],mom2[0],mom2[1],mom2[2],val.real(),val.imag() };
    fwrite(dbuf, sizeof(double), 9 ,fp);
  }
}


void _ContractedWallSinkBilinearSpecMomentum_helper<SpinColorFlavorMatrix>::write(std::ostream &fp, const Rcomplex &val, const int &idx1, const int &idx2, const Float &p2, 
										  const std::vector<Float> &mom1, const std::vector<Float> &mom2, const int &t){
  std::pair<int,int> gamma_sigma1 = _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unmap(idx1);
  std::pair<int,int> gamma_sigma2 = _PropagatorBilinear_helper<SpinColorFlavorMatrix>::unmap(idx2);

  //format is gamma1 sigma1 gamma2 sigma2 p^2 p1x p1y p1z p2x p2y p2z t val.re val.im
  if(!UniqueID()) 
    fp << gamma_sigma1.first << " " << gamma_sigma1.second << " " 
       << gamma_sigma2.first << " " << gamma_sigma2.second << " "
       << p2 << " " << mom1[0] << " " << mom1[1] << " " << mom1[2] << " " << mom2[0] << " " << mom2[1] << " " << mom2[2] << " " << t << " " << val.real() << " " << val.imag() << std::endl;
}


void _ContractedWallSinkBilinearSpecMomentum_helper<WilsonMatrix>::write(std::ostream &fp, const Rcomplex &val, const int &idx1, const int &idx2, const Float &p2, 
									 const std::vector<Float> &mom1, const std::vector<Float> &mom2, const int &t){
  //format is gamma1 gamma2 p^2 p1x p1y p1z p2x p2y p2z t val.re val.im
  if(!UniqueID())
    fp << idx1 << " " << idx2 << " " << p2 << " " << mom1[0] << " " << mom1[1] << " " << mom1[2] << " " << mom2[0] << " " << mom2[1] << " " << mom2[2] << " " << t << " " << val.real() << " " << val.imag() << std::endl;
}


void doFlipMomentum<WilsonMatrix>::flip(WilsonMatrix &mat){
  ERR.General("doFlipMomentum<WilsonMatrix>","flip","Only valid for G-parity BCs\n");
}
void doFlipMomentum<SpinColorFlavorMatrix>::flip(SpinColorFlavorMatrix &mat){
  mat.flipSourceMomentum();
}

CPS_END_NAMESPACE