C++ Lattice::GaugeField方法代码示例

// multiply exp( - i Q ) in `mu' direction. 
void twist_links(Lattice &lat, const Float Q, const int mu)
{
  // multiply su3 links by u1 links
  // assumes both are in canonical order!
  // link =  exp(-i A_mu Q) where Q is the 
  // quark electric charge or twist angle
  int x[4];
  Complex cc(cos(Q), -sin(Q)); 
  
  Matrix temp;
  for(x[3]=0; x[3]<GJP.TnodeSites();x[3]++){
    for(x[2]=0; x[2]<GJP.ZnodeSites();x[2]++){
      for(x[1]=0; x[1]<GJP.YnodeSites();x[1]++){
	for(x[0]=0; x[0]<GJP.XnodeSites();x[0]++){
	  int offset_x = lat.GsiteOffset(x);
	  // the link matrix at current site, in direction mu:
	  Matrix *su3_link = lat.GaugeField() + offset_x + mu;
	  //cTimesVec((float*)&temp, *(float*)phase, *((float*)phase+1), 
	  cTimesVec((IFloat*)&temp, (IFloat)cc.real(), (IFloat)cc.imag(), 
		    (IFloat*)su3_link, 18);
	  moveMem((IFloat*)su3_link, (IFloat*)&temp, 18);
	}
      }
    }
  }
}

void pt_init(Lattice &Lat){
  PTArg pt_arg;
//Size of local volume in all four directions
  pt_arg.size[0] = GJP.XnodeSites();
  pt_arg.size[1] = GJP.YnodeSites();
  pt_arg.size[2] = GJP.ZnodeSites();
  pt_arg.size[3] = GJP.TnodeSites();
  //-------------------------------------------------
  //Added by Michael C.
  for(int i = 0; i < 4;i++)
    if(GJP.Nodes(i) == 1)
      pt_arg.local[i] = 1;
    else
      pt_arg.local[i] = 0;
  //-------------------------------------------------
  int temp = 0;
  for(int i = 0;i<4;i++){
    temp += GJP.NodeSites(i)*GJP.NodeCoor(i);
  }
  if (temp%2==0)
    pt_arg.evenodd = PT_EVEN;
  else
    pt_arg.evenodd = PT_ODD;
  pt_arg.gauge_field_addr = (IFloat *)Lat.GaugeField();

  StrOrdType str_ord = Lat.StrOrd();
  printf("str_ord=%d\n",str_ord);

  switch(str_ord){
    case CANONICAL:
      pt_arg.g_str_ord = PT_XYZT;
      pt_arg.v_str_ord = PT_XYZT;
      pt_arg.v_str_ord_cb = PT_TXYZ;
      pt_arg.g_conj = 0;
      break;
    case WILSON:
      pt_arg.g_str_ord = PT_XYZT_CB_O;
      pt_arg.v_str_ord = PT_XYZT_CB_O;
      pt_arg.v_str_ord_cb = PT_TXYZ;
      pt_arg.g_conj = 0;
      break;
    case STAG:
      pt_arg.g_str_ord = PT_XYZT;
      pt_arg.v_str_ord = PT_XYZT;
      pt_arg.v_str_ord_cb = PT_TXYZ;
      pt_arg.g_conj = 1;
      break;
    default:
      break;
  }
  pt_arg.prec = sizeof(IFloat);
  StaticPT.init(&pt_arg);
}

void PlaquetteMonster::Initialize(const Lattice & lat)
{

  Copy( lat, *lats[0]);

  int s[4];
  int t[4];
  for ( int m=1; m<len; ++m) {
    Copy( *lats[m-1], *lats[m]);
    for ( s[0]=0, t[0]=0; s[0]<sites[0]; ++s[0], ++t[0] )
    for ( s[1]=0, t[1]=0; s[1]<sites[1]; ++s[1], ++t[1] )
    for ( s[2]=0, t[2]=0; s[2]<sites[2]; ++s[2], ++t[2] )
    for ( s[3]=0, t[3]=0; s[3]<sites[3]; ++s[3], ++t[3] )
    for (int mu=0; mu<4; ++mu) {
      t[mu]+=m;
      lats[m]->GaugeField(mu,s) *= lat.GaugeField(mu,t);
      t[mu]-=m;
    }
  }
}

// ----------------------------------------------------------------
// twisted_bc: set twisted boundary condition.
//
// The gauge field must be in CANONICAL order.
//
// add: if true then multiply exp(i*angle), otherwise multiply
// exp(-i*angle). This parameter can be used to add/remove the bc.
// ----------------------------------------------------------------
void twisted_bc(Lattice &lat, const double mom[4], bool add)
{
    Matrix *gauge = (Matrix *)(lat.GaugeField());

    const double PI = 3.1415926535897932384626433832795028842;
    double sign = add ? 1 : -1;

    for(int mu = 0; mu < 4; ++mu) {
        if(mom[mu] == 0) continue;
        // if(GJP.NodeCoor(mu) + 1 != GJP.Nodes(mu)) continue;

        double t = 2.0 * PI / GJP.Sites(mu) * mom[mu];
        const Rcomplex cval(cos(t), sign * sin(t));

        int low[4] = { 0, 0, 0, 0 };
        int high[4] = { GJP.XnodeSites(), GJP.YnodeSites(),
                        GJP.ZnodeSites(), GJP.TnodeSites() };
        // low[mu] = high[mu] - 1;

        int hl[4] = { high[0] - low[0], high[1] - low[1],
                      high[2] - low[2], high[3] - low[3] };

        const int hl_sites = hl[0] * hl[1] * hl[2] * hl[3];

#pragma omp parallel for
        for(int i = 0; i < hl_sites; ++i) {
            int x[4];
            compute_coord(x, hl, low, i);

            int off = mu + 4 * compute_id(x, high);
            gauge[off] *= cval;
        }
    }

    Fbfm *fbfm = dynamic_cast<Fbfm *>(&lat);
    if(fbfm != NULL) {
        fbfm->ImportGauge();
    }
}

CPS_START_NAMESPACE

#define PROFILE
void WriteLatticeParallel::write(Lattice & lat, const QioArg & wt_arg)
{
  const char * fname = "write()";
  VRB.Func(cname,fname);

  char loginfo[100];
  sprintf(loginfo,"Unload %s",wt_arg.FileName);
  startLogging(loginfo);

#ifdef PROFILE
  struct timeval start,end;
  gettimeofday(&start,NULL);
#endif
  //sync();

  // init
  int error = 0;
  io_good = false;

  //  const char * filename = wt_arg.FileName;
  recon_row_3 = wt_arg.ReconRow3;
  FP_FORMAT dataFormat = wt_arg.FileFpFormat;
  fpconv.setFileFormat(dataFormat);
  if(fpconv.fileFormat == FP_UNKNOWN) {
    ERR.General(cname,fname,"Output Floating Point format UNKNOWN\n");
  }

  // calc Plaq and LinkTrace
  const int size_matrices(wt_arg.VolNodeSites()*4);
  Matrix * lpoint = lat.GaugeField();
  VRB.Flow(cname,fname, "Writing Gauge Field at Lattice::GaugeField() = %p\n", lpoint);

  Float plaq = lat.SumReTrPlaq()/(18*wt_arg.VolSites()) ;
  Float ltrace(0.0);
  if(wt_arg.Scoor() == 0) {
    for(int i=0;i<size_matrices;i++){
      ltrace += (lpoint+i)->ReTr();
    }
    ltrace = globalSumFloat(ltrace) / (4*3*wt_arg.VolSites());
  }
  else
    globalSumFloat(0.0);  // everyone has to participate in global ops

  log();
  
  // write lattice data, in Parallel or Serial manner
  // determined by the template parameter "IoStyle" of this class
  int data_per_site = recon_row_3 ? 4*12 : 4*18;
  const int chars_per_site = data_per_site * fpconv.fileFpSize();

  unsigned int csum = 0;

#if TARGET != QCDOC   // when not on QCDOC(like on LINUX), use serial IO mode
  setSerial();
#endif
  
  fstream output;

  if (isRoot()){
    output.open(wt_arg.FileName,ios::out);
	output.close();
  }
  Float temp=0.;
  glb_sum(&temp);
  //sync();
  if(parIO()) {
    // all open file, start writing
    output.open(wt_arg.FileName);
    if(!output.good())    {
            ERR.General(cname,fname, "Could not open file: [%s] for output.\n",wt_arg.FileName);
      //      VRB.Flow(cname,fname,"USER: maybe you should kill the process\n");
      
      printf("Node %d:Could not open file: [%s] for output.\n",UniqueID(),wt_arg.FileName);
      error = 1;
    }
  }
  else {
    // only node 0 open file, start writing
    if(isRoot()) {
      FILE *fp = Fopen(wt_arg.FileName,"w");
      Fclose(fp);
      output.open(wt_arg.FileName);
      if(!output.good())    {
	//	VRB.Flow(cname,fname, "Could not open file: [%s] for output.\n",wt_arg.FileName);
	//	VRB.Flow(cname,fname,"USER: maybe you should kill the process\n");
	error = 1;
      }
    }
  }
  if (error)
    printf("Node %d: says opening %s failed\n",UniqueID(),wt_arg.FileName);
//  if(synchronize(error) > 0)  
//    ERR.FileW(cname,fname,wt_arg.FileName);
   error=0;

  // write header
  if(isRoot()){
//.........这里部分代码省略.........

static int cps_qdp_dwf_invert(
Lattice &Lat, 
CPS_NAMESPACE::Float mass, 
CPS_NAMESPACE::Float *sol_cps, 
CPS_NAMESPACE::Float *src_cps,
double residual,
int max_iter
){
  char *fname = "cps_qdp_dwf_invert()";

  if (GJP.Snodes() != 1)
    ERR.General("",fname,"Snodes()(%d)!=1",GJP.Snodes());

  cps_qdp_init(GJP.argc_p(), GJP.argv_p());
  if (!qdp_initted)
    ERR.General("",fname,"QDP not initialized!(%d)",0);

  int iter=0;   

  int Ls = GJP.SnodeSites();
  double M5 = GJP.DwfHeight();
//  int Nd = 4;

  Printf("src[0]=%g\n",*src_cps);
  int lx = QDP::Layout::subgridLattSize()[0];
  int ly = QDP::Layout::subgridLattSize()[1];
  int lz = QDP::Layout::subgridLattSize()[2];
  int lt = QDP::Layout::subgridLattSize()[3];
  multi1d<LatticeFermion> src_qdp(Ls);

  /********************************************************
   * Setup DWF operator
   ********************************************************
   */
  omp_set_num_threads(64);
  bfmarg dwfa;
#ifdef UNIFORM_SEED_TESTING
  majorityVote  dwf;
#else
//  bfm_qdp<double>  dwf;
  bfm_qdp<float>  dwf;
#endif

  dwfa.node_latt[0]  = lx;
  dwfa.node_latt[1]  = ly;
  dwfa.node_latt[2]  = lz;
  dwfa.node_latt[3]  = lt;
  dwfa.verbose=1;
  dwfa.reproduce=0;
  bfmarg::Threads(64);
  bfmarg::Reproduce(0);
  bfmarg::ReproduceChecksum(0);
  bfmarg::ReproduceMasterCheck(0);
  bfmarg::Verbose(1);

  static int CGcount = 0;
  CGcount++;
  int test_freq = GJP.CGreprodFreq();
  if (test_freq && (CGcount % test_freq == 0)) {
  dwfa.reproduce=1;
  bfmarg::Reproduce(1);
  }

  multi1d<int> procs = QDP::Layout::logicalSize();
  Printf("%d dim machine\n\t", procs.size());
  for(int mu=0;mu<4;mu++){
    Printf("%d ", procs[mu]);
    if ( procs[mu]>1 ) {
      dwfa.local_comm[mu] = 0;
    } else { 
      dwfa.local_comm[mu] = 1;
    }
  }
  Printf("\nLocal comm = ");
  for(int mu=0;mu<4;mu++){
    Printf("%d ", dwfa.local_comm[mu]);
  }
  Printf("\n");
  
  multi1d<int> ncoor = QDP::Layout::nodeCoord();

//  Real M5(1.8);
//  Real mq(0.1);

  CPS_NAMESPACE::Float *gauge = (CPS_NAMESPACE::Float*) Lat.GaugeField();
  dwfa.precon_5d = 1;
  dwfa.Ls   = Ls;
  dwfa.solver = DWF;
  dwfa.M5   = toDouble(M5);
  dwfa.mass = toDouble(mass);
  dwfa.Csw  = 0.0;
  dwfa.max_iter = max_iter;
  dwfa.residual = residual;
//OK
for(int i = 0;i<1;i++){
  Printf("Initialising bfm operator\n");
  dwf.init(dwfa);
  Printf("Initialising bfm operator done\n");

//OK
//.........这里部分代码省略.........

void LatticeContainer::Set(Lattice &lat){
	if (str_ord != lat.StrOrd())
	ERR.General(cname,"Set()","Storage ordering of LatticeContainer(%d) doesn't agree with lattice ordering(%d)", str_ord,lat.StrOrd());
	lat.GaugeField(gauge_p);
}

//.........这里部分代码省略.........
  //-------------------------------------------------------------------------------

  int tt,xx,yy,zz,temp_ind;
  int slice_ind[3];
  
  //slice_ind[3] stores the 3 directions on the 'dir' slice with indices increasing 
  temp_ind = 0;
  for (int i=0; i<4; i++){
    if (i!=dir) {
      slice_ind[temp_ind] = i;
      temp_ind++;
    }
  }
  
  VRB.Result(cname,fname,"dir == %d \n",dir);
  VRB.Result(cname,fname,"slice index == %d %d %d\n",slice_ind[0],slice_ind[1],slice_ind[2]);
  
  // the dummy node_sites for each dummy dirction
  int NN[4];
  NN[0] = node_sites[slice_ind[0]];
  NN[1] = node_sites[slice_ind[1]];
  NN[2] = node_sites[slice_ind[2]];
  NN[3] = node_sites[dir];
  
  int s[4];
  for (int i=0; i<3; i++)
    s[i] = slice_ind[i];
  s[3] = dir;

  //---------------------------------------------------------------------------------
  //copy the old lattice config to matrix array L, transform L and then copy back
  //---------------------------------------------------------------------------------
  
  lat.CopyGaugeField(L);  
  int x[4];
  
  // xx yy zz tt are dummy position vector, as tt represents the gfixing direction  
  for (x[3]=0; x[3]<node_sites[3]; x[3]++)
    for (x[2]=0; x[2]<node_sites[2]; x[2]++)
      for (x[1]=0; x[1]<node_sites[1]; x[1]++)
	for (x[0]=0; x[0]<node_sites[0]; x[0]++)
	  {
	    xx = x[slice_ind[0]];
	    yy = x[slice_ind[1]];
	    zz = x[slice_ind[2]];
	    tt = x[dir];
	
	    Matrix g = Gp[tt][(zz*NN[1]+yy)*NN[0]+xx];
	    Matrix D; 
	    
	    Matrix gg ;
	    Matrix transmit;
	    
	    //----------------- T Direction ----------------------------
	    if (tt+1<NN[3]) gg = Gp[tt+1][(zz*NN[1]+yy)*NN[0]+xx];
	    else { 
	      transmit = Gp[0][(zz*NN[1]+yy)*NN[0]+xx];
	      getPlusData((IFloat *)&gg, (IFloat *)&transmit,
			  sizeof(Matrix)/sizeof(IFloat), dir) ;
	    }
	    D.Dagger(gg);
	    L[IND(x[0],x[1],x[2],x[3],dir)] = g*L[IND(x[0],x[1],x[2],x[3],dir)]*D;
	    
	    //----------------- Z Direction ----------------------------
	    if (zz+1<NN[2]) gg = Gp[tt][((zz+1)*NN[1]+yy)*NN[0]+xx]; 
	    else {