C++ Matrix::Dagger方法代码示例

void pt_mat(int N, IFloat **fout, IFloat **fin, const int *dir){

    IFloat *u, *v;
    int s[4], spm[4];
    Matrix ud;
    Matrix *vin[N],*vout[N];

    for(int d=0; d<N; d++){
	vin[d] = (Matrix *)fin[d];
	vout[d] = (Matrix *)fout[d];

	int mu = dir[d];
	bool backwards = mu%2 ? true : false;
        mu /= 2;

	for(spm[3]=s[3]=0; s[3]<GJP.TnodeSites(); spm[3]++, s[3]++)
	    for(spm[2]=s[2]=0; s[2]<GJP.ZnodeSites(); spm[2]++, s[2]++)
		for(spm[1]=s[1]=0; s[1]<GJP.YnodeSites(); spm[1]++, s[1]++)
		    for(spm[0]=s[0]=0; s[0]<GJP.XnodeSites(); spm[0]++, s[0]++) {
			if(backwards){
			    spm[mu] = (s[mu]-1+GJP.NodeSites(mu))%GJP.NodeSites(mu);
			    u = (IFloat*)(Lat->GaugeField()+Lat->GsiteOffset(spm)+mu);
			    v = (IFloat*)&vin[d][Lat->FsiteOffset(spm)];
			    mDotMEqual((IFloat*)&vout[d][Lat->FsiteOffset(s)],
				       u, v);
			}else{
			    spm[mu] = (s[mu]+1)%GJP.NodeSites(mu);
			    ud.Dagger(*(Lat->GaugeField()+Lat->GsiteOffset(s)+mu));
			    u = (IFloat*)&ud;
			    v = (IFloat*)&vin[d][Lat->FsiteOffset(spm)];
			    // STAG stores hermitian conjugate links
			    mDotMEqual((IFloat*)&vout[d][Lat->FsiteOffset(s)],
				      u, v);			    
			}
			spm[mu] = s[mu];

		    }
			
    }
    
}

CPS_START_NAMESPACE

//----------------------------------------------------------
//
// alg_tcharge.C
// 
// measures a simple defintion of the
// topological charge
//  
//  Using Clover Leaf   c.f   hep-lat/010610  Eq (6) -- (10)
//-----------------------------------------------------------


/*!
  takes the imaginary part of a matrix
*/
void ZeroReal(Matrix& m)
{
  Matrix temp;
  temp.Dagger( m );
  m-=temp;
  m*=0.5;
}

//.........这里部分代码省略.........
  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 {
	      transmit = Gp[tt][((zz+1)%NN[2]*NN[1]+yy)*NN[0]+xx]; 
	      getPlusData((IFloat *)&gg, (IFloat *)&transmit,
			  sizeof(Matrix)/sizeof(IFloat), slice_ind[2]) ;
	    }
	    D.Dagger(gg);
	    L[IND(x[0],x[1],x[2],x[3],slice_ind[2])] = g*L[IND(x[0],x[1],x[2],x[3],slice_ind[2])]*D;
	    
	    //----------------- Y Direction ----------------------------
	    if (yy+1<NN[1]) gg = Gp[tt][(zz*NN[1]+(yy+1))*NN[0]+xx];
	    else {
	      transmit = Gp[tt][(zz*NN[1]+(yy+1)%NN[1])*NN[0]+xx];
	      getPlusData((IFloat *)&gg, (IFloat *)&transmit,
			  sizeof(Matrix)/sizeof(IFloat), slice_ind[1]) ;
	    }
	    D.Dagger(gg);
	    L[IND(x[0],x[1],x[2],x[3],slice_ind[1])] = g*L[IND(x[0],x[1],x[2],x[3],slice_ind[1])]*D;
	    
	    //----------------- X Direction ----------------------------
	    if (xx+1<NN[0]) gg = Gp[tt][(zz*NN[1]+yy)*NN[0]+(xx+1)];
	    else {
	      transmit = Gp[tt][(zz*NN[1]+yy)*NN[0]+(xx+1)%NN[0]];
	      getPlusData((IFloat *)&gg, (IFloat *)&transmit,
			  sizeof(Matrix)/sizeof(IFloat), slice_ind[0]) ;
	    }
	    D.Dagger(gg);
	    L[IND(x[0],x[1],x[2],x[3],slice_ind[0])] = g*L[IND(x[0],x[1],x[2],x[3],slice_ind[0])]*D;
	  }
  
  
  lat.GaugeField(L);
  
  //--------------------------------Free L and Gp -------------------------------------
  sfree(L);
  
  for (ii=0; ii<npln; ii++)
    sfree(Gp[ii]);
  
  sfree(Gp);
  sfree(plns);

}

//.........这里部分代码省略.........
      aQ=seqQ[s.Index()] ;
      // Apply the gamma matrices on the sequential  propagator
      // We multiply the source indices of the sequential propagator
      // which are the anti-quark
      for(int mu(0);mu<G.N();mu++)
	aQ.gr(G[mu]) ; // Multiply by gamma_mu 

      // Loop over all derivative terms on the Quark
      for(D.Start();D.NotEnd();D.NextQuark()) 
	{
	  // Calculate the endpoints of the shift
	  D.CalcEndPoints(aq_vec,q_vec,s);
	  // The sink indices of the quark propagator are the quark
	  Q=Quark.GetMatrix(q_vec,buff) ;
	  
	 
	  D.CalcListOfLinks(dirs);

#ifdef DEBUG
printf("\n%s:: dirs ",cname);
for(int i(0);i<D.N();i++) printf("%i ",dirs[i]);
printf(" aq_vec: ");
for(int i(0);i<4;i++) printf("%i ",aq_vec[i]);
printf(" q_vec: ") ; 
for(int i(0);i<4;i++) printf("%i ",q_vec[i]);
#endif

	  // Start from the anti-quark and move along the path described
	  // by the list of dirs multiplying all the encountered links
	  // at the end of the day you should be at the position of the quark
          if(D.N()>1)  
	    {
	      U.ZeroMatrix();
	      lat.PathOrdProdPlus(U,aq_vec,dirs,D.N()) ;
	      //lat.PathOrdProd(U,aq_vec,dirs,D.N()) ;
	    }
	  else //if D.N()==1 PathOrdPlus does not work
	    {
	      int abs_dir = dirs[0]&3; //abs_dir is {0,1,2,3}
	      if(abs_dir == dirs[0]) // dir is forward
		U = *(lat.GetBufferedLink(aq_vec, abs_dir));
	      else // dir is backward. q_vec is now the site where the link is
		U.Dagger(*(lat.GetBufferedLink(q_vec, abs_dir)));
	    }

#ifdef DEBUG
Complex cboo(U.Tr()) ; 
printf("  TrU= (%g %g)\n",cboo.real(),cboo.imag());       
#endif

	  //Now that the gauge connection is constructed multiply it on
	  //the quark i.e the source indices of the anti-quark propagator
	  aQU.UMultSource(U,aQ) ;

	  Terms(D.DTermIndx(),s.Index()) = Trace(aQU,Q) ;

#ifdef DEBUG
printf("   (anti-Q, Q) aQU = (%6.3f %6.3f), (%6.3f %6.3f) (%6.3f %6.3f)\n",
       aQ.Trace().real(),aQ.Trace().imag(),
       Q.Trace().real(),Q.Trace().imag(),
       aQU.Trace().real(),aQU.Trace().imag());
printf("%s:: DTermIndx(): %i site: %i Value: ( %6.3f %6.3f ) \n",
       cname,D.DTermIndx(),s.Index(),
       Terms(D.DTermIndx(),s.Index()).real(), 
       Terms(D.DTermIndx(),s.Index()).imag());
#endif


	}
    }

  Complex tt ;
  // Loop over sites and construct the operator
  for(s.Begin();s.End();s.nextSite()) 
    {
      int t(s.physT()) ;
      for(D.Start();D.NotEnd();D.Next()) // Loop over all derivative terms
	{
	  //Calculate the endpoints of the shift
	  //aq_vec : anti quark location
	  // q_vec :      quark location
	  D.CalcEndPoints(aq_vec,q_vec,s);
	  Complex cc(Terms.GetTerm(D.DTermIndx(),aq_vec,tt)) ;

	  //Finaly the momentum factors
	  cc *= mom.Fact(s) ; // The momentum factor exp(-ipx)	
	  cc *= D.Fact()    ; /* The sign and the normalization factor 
				 for the current derivative term         */

	  tmp[t] += cc ;

#ifdef DEBUG
printf("%s:: DTermIndx(): %i site: %i Fact: %9.6f Value: ( %6.3f %6.3f ) \n",
       cname, D.DTermIndx(),s.Index(),D.Fact(),cc.real(),cc.imag());
#endif
	
	}
    }
  delete [] dirs ;
}