mirror of
				https://github.com/paboyle/Grid.git
				synced 2025-11-03 21:44:33 +00:00 
			
		
		
		
	Introduce view objects that can safely be copied to GPU for access
This commit is contained in:
		@@ -109,12 +109,14 @@ int main(int argc, char ** argv) {
 | 
				
			|||||||
	  StencilEntry *SE;
 | 
						  StencilEntry *SE;
 | 
				
			||||||
	  SE = myStencil.GetEntry(permute_type,0,i);
 | 
						  SE = myStencil.GetEntry(permute_type,0,i);
 | 
				
			||||||
	  
 | 
						  
 | 
				
			||||||
 | 
						  auto check = Check.View();
 | 
				
			||||||
 | 
						  auto foo   = Foo.View();
 | 
				
			||||||
	  if ( SE->_is_local && SE->_permute )
 | 
						  if ( SE->_is_local && SE->_permute )
 | 
				
			||||||
	    permute(Check[i],Foo[SE->_offset],permute_type);
 | 
						    permute(check[i],foo[SE->_offset],permute_type);
 | 
				
			||||||
	  else if (SE->_is_local)
 | 
						  else if (SE->_is_local)
 | 
				
			||||||
	    Check[i] = Foo[SE->_offset];
 | 
						    check[i] = foo[SE->_offset];
 | 
				
			||||||
	  else { 
 | 
						  else { 
 | 
				
			||||||
	    Check[i] = myStencil.CommBuf()[SE->_offset];
 | 
						    check[i] = myStencil.CommBuf()[SE->_offset];
 | 
				
			||||||
	    //	    std::cout << " receive "<<i<<" " << Check[i]<<std::endl;
 | 
						    //	    std::cout << " receive "<<i<<" " << Check[i]<<std::endl;
 | 
				
			||||||
	    //	    std::cout << " Foo     "<<i<<" " <<   Foo[i]<<std::endl;
 | 
						    //	    std::cout << " Foo     "<<i<<" " <<   Foo[i]<<std::endl;
 | 
				
			||||||
	  }
 | 
						  }
 | 
				
			||||||
@@ -149,8 +151,10 @@ int main(int argc, char ** argv) {
 | 
				
			|||||||
	}}}}
 | 
						}}}}
 | 
				
			||||||
 | 
					
 | 
				
			||||||
	if (nrm > 1.0e-4) {
 | 
						if (nrm > 1.0e-4) {
 | 
				
			||||||
	  for(int i=0;i<Check.size();i++){
 | 
						  auto check = Check.View();
 | 
				
			||||||
	    std::cout << i<<" Check.odata "<<Check[i]<< "\n"<<i<<" Bar.odata "<<Bar[i]<<std::endl;
 | 
						  auto bar   = Bar.View();
 | 
				
			||||||
 | 
						  for(int i=0;i<check.size();i++){
 | 
				
			||||||
 | 
						    std::cout << i<<" Check "<<check[i]<< "\n"<<i<<" Bar "<<bar[i]<<std::endl;
 | 
				
			||||||
	  }
 | 
						  }
 | 
				
			||||||
	}
 | 
						}
 | 
				
			||||||
	if (nrm > 1.0e-4) exit(-1);
 | 
						if (nrm > 1.0e-4) exit(-1);
 | 
				
			||||||
@@ -206,12 +210,14 @@ int main(int argc, char ** argv) {
 | 
				
			|||||||
	  SE = EStencil.GetEntry(permute_type,0,i);
 | 
						  SE = EStencil.GetEntry(permute_type,0,i);
 | 
				
			||||||
	  //	  std::cout << "Even source "<< i<<" -> " <<SE->_offset << " "<< SE->_is_local<<std::endl;
 | 
						  //	  std::cout << "Even source "<< i<<" -> " <<SE->_offset << " "<< SE->_is_local<<std::endl;
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 | 
						  auto ocheck = OCheck.View();
 | 
				
			||||||
 | 
						  auto efoo   = EFoo.View();
 | 
				
			||||||
	  if ( SE->_is_local && SE->_permute )
 | 
						  if ( SE->_is_local && SE->_permute )
 | 
				
			||||||
	    permute(OCheck[i],EFoo[SE->_offset],permute_type);
 | 
						    permute(ocheck[i],efoo[SE->_offset],permute_type);
 | 
				
			||||||
	  else if (SE->_is_local)
 | 
						  else if (SE->_is_local)
 | 
				
			||||||
	    OCheck[i] = EFoo[SE->_offset];
 | 
						    ocheck[i] = efoo[SE->_offset];
 | 
				
			||||||
	  else
 | 
						  else
 | 
				
			||||||
	    OCheck[i] = EStencil.CommBuf()[SE->_offset];
 | 
						    ocheck[i] = EStencil.CommBuf()[SE->_offset];
 | 
				
			||||||
	}
 | 
						}
 | 
				
			||||||
	OStencil.HaloExchange(OFoo,compress);
 | 
						OStencil.HaloExchange(OFoo,compress);
 | 
				
			||||||
	for(int i=0;i<ECheck.Grid()->oSites();i++){
 | 
						for(int i=0;i<ECheck.Grid()->oSites();i++){
 | 
				
			||||||
@@ -220,12 +226,14 @@ int main(int argc, char ** argv) {
 | 
				
			|||||||
	  SE = OStencil.GetEntry(permute_type,0,i);
 | 
						  SE = OStencil.GetEntry(permute_type,0,i);
 | 
				
			||||||
	  //	  std::cout << "ODD source "<< i<<" -> " <<SE->_offset << " "<< SE->_is_local<<std::endl;
 | 
						  //	  std::cout << "ODD source "<< i<<" -> " <<SE->_offset << " "<< SE->_is_local<<std::endl;
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 | 
						  auto echeck = ECheck.View();
 | 
				
			||||||
 | 
						  auto ofoo   = OFoo.View();
 | 
				
			||||||
	  if ( SE->_is_local && SE->_permute )
 | 
						  if ( SE->_is_local && SE->_permute )
 | 
				
			||||||
	    permute(ECheck[i],OFoo[SE->_offset],permute_type);
 | 
						    permute(echeck[i],ofoo[SE->_offset],permute_type);
 | 
				
			||||||
	  else if (SE->_is_local)
 | 
						  else if (SE->_is_local)
 | 
				
			||||||
	    ECheck[i] = OFoo[SE->_offset];
 | 
						    echeck[i] = ofoo[SE->_offset];
 | 
				
			||||||
	  else
 | 
						  else
 | 
				
			||||||
	    ECheck[i] = OStencil.CommBuf()[SE->_offset];
 | 
						    echeck[i] = OStencil.CommBuf()[SE->_offset];
 | 
				
			||||||
	}
 | 
						}
 | 
				
			||||||
 | 
					
 | 
				
			||||||
	setCheckerboard(Check,ECheck);
 | 
						setCheckerboard(Check,ECheck);
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -30,7 +30,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 | 
				
			|||||||
#include <Grid/Grid.h>
 | 
					#include <Grid/Grid.h>
 | 
				
			||||||
 | 
					
 | 
				
			||||||
using namespace Grid;
 | 
					using namespace Grid;
 | 
				
			||||||
 ;
 | 
					 | 
				
			||||||
 | 
					
 | 
				
			||||||
#define POWER10
 | 
					#define POWER10
 | 
				
			||||||
 | 
					
 | 
				
			||||||
@@ -84,10 +83,11 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
  pickCheckerboard(Odd,Uo,U);
 | 
					  pickCheckerboard(Odd,Uo,U);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
  //  std::cout<<GridLogMessage << U<<std::endl;
 | 
					  //  std::cout<<GridLogMessage << U<<std::endl;
 | 
				
			||||||
 | 
					  std::cout<<GridLogMessage<< U <<std::endl;
 | 
				
			||||||
 | 
					  std::cout<<GridLogMessage << "U " <<norm2(U)<<std::endl;
 | 
				
			||||||
  std::cout<<GridLogMessage << "Ue " <<norm2(Ue)<<std::endl;
 | 
					  std::cout<<GridLogMessage << "Ue " <<norm2(Ue)<<std::endl;
 | 
				
			||||||
  std::cout<<GridLogMessage << "Uo " <<norm2(Uo)<<std::endl;
 | 
					  std::cout<<GridLogMessage << "Uo " <<norm2(Uo)<<std::endl;
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 | 
					 | 
				
			||||||
  TComplex cm;
 | 
					  TComplex cm;
 | 
				
			||||||
  TComplex cmeo;
 | 
					  TComplex cmeo;
 | 
				
			||||||
  for(int dir=0;dir<Nd;dir++){
 | 
					  for(int dir=0;dir<Nd;dir++){
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -547,7 +547,7 @@ int main(int argc, char **argv) {
 | 
				
			|||||||
                   coor[1]++) {
 | 
					                   coor[1]++) {
 | 
				
			||||||
                for (coor[0] = 0; coor[0] < latt_size[0] / mpi_layout[0];
 | 
					                for (coor[0] = 0; coor[0] < latt_size[0] / mpi_layout[0];
 | 
				
			||||||
                     coor[0]++) {
 | 
					                     coor[0]++) {
 | 
				
			||||||
                  std::complex<Grid::Real> diff;
 | 
					                  Complex diff;
 | 
				
			||||||
 | 
					
 | 
				
			||||||
                  std::vector<int> shiftcoor = coor;
 | 
					                  std::vector<int> shiftcoor = coor;
 | 
				
			||||||
                  shiftcoor[dir] = (shiftcoor[dir] + shift + latt_size[dir]) %
 | 
					                  shiftcoor[dir] = (shiftcoor[dir] + shift + latt_size[dir]) %
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -89,9 +89,11 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
  // replicate across fifth dimension
 | 
					  // replicate across fifth dimension
 | 
				
			||||||
  ////////////////////////////////////
 | 
					  ////////////////////////////////////
 | 
				
			||||||
  LatticeGaugeField Umu5d(FGrid); 
 | 
					  LatticeGaugeField Umu5d(FGrid); 
 | 
				
			||||||
 | 
					  auto umu5d = Umu5d.View();
 | 
				
			||||||
 | 
					  auto umu   = Umu.View();
 | 
				
			||||||
  for(int ss=0;ss<Umu.Grid()->oSites();ss++){
 | 
					  for(int ss=0;ss<Umu.Grid()->oSites();ss++){
 | 
				
			||||||
    for(int s=0;s<Ls;s++){
 | 
					    for(int s=0;s<Ls;s++){
 | 
				
			||||||
      Umu5d[Ls*ss+s] = Umu[ss];
 | 
					      umu5d[Ls*ss+s] = umu[ss];
 | 
				
			||||||
    }
 | 
					    }
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -98,15 +98,18 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
					    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    // fourth order exponential approx
 | 
					    // fourth order exponential approx
 | 
				
			||||||
    thread_loop( (auto i=mom.begin();i<mom.end();i++),{
 | 
					    auto mom_v = mom.View();
 | 
				
			||||||
      Uprime[i](mu) =
 | 
					    auto Uprime_v = Uprime.View();
 | 
				
			||||||
	  U[i](mu)
 | 
					    auto U_v = U.View();
 | 
				
			||||||
	+ mom[i](mu)*U[i](mu)*dt 
 | 
					    thread_loop( (auto i=mom_v.begin();i<mom_v.end();i++),{
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
 | 
					      Uprime_v[i](mu) =
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
 | 
						  U_v[i](mu)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
 | 
						+ mom_v[i](mu)*U_v[i](mu)*dt 
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
				
			||||||
	;
 | 
						;
 | 
				
			||||||
    });
 | 
					    });
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -112,29 +112,35 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
    // fourth order exponential approx
 | 
					    // fourth order exponential approx
 | 
				
			||||||
 | 
					
 | 
				
			||||||
#if 0
 | 
					#if 0
 | 
				
			||||||
    thread_loop( (auto i=mom.begin();i<mom.end();i++),{
 | 
					    auto mom_v    = mom.View();
 | 
				
			||||||
      Uprime[i](mu) =
 | 
					    auto U_v      = U.View();
 | 
				
			||||||
	  U[i](mu)
 | 
					    auto Uprime_v = Uprime.View();
 | 
				
			||||||
	+ mom[i](mu)*U[i](mu)*dt 
 | 
					
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
 | 
					    thread_loop( (auto i=mom_v.begin();i<mom_v.end();i++),{
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
 | 
					      Uprime_v[i](mu) =	  U_v[i](mu)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
 | 
						+ mom_v[i](mu)*U_v[i](mu)*dt 
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
				
			||||||
	;
 | 
						;
 | 
				
			||||||
      });
 | 
					    });
 | 
				
			||||||
#else
 | 
					#else
 | 
				
			||||||
    accelerator_loop( i, mom, {
 | 
					    auto mom_v    = mom.View();
 | 
				
			||||||
      Uprime[i](mu) =
 | 
					    auto U_v      = U.View();
 | 
				
			||||||
	  U[i](mu)
 | 
					    auto Uprime_v = Uprime.View();
 | 
				
			||||||
	+ mom[i](mu)*U[i](mu)*dt 
 | 
					
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
 | 
					    accelerator_loop( i,mom_v,{
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
 | 
					      Uprime_v[i](mu) =	  U_v[i](mu)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
 | 
						+ mom_v[i](mu)*U_v[i](mu)*dt 
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
				
			||||||
	;
 | 
						;
 | 
				
			||||||
      });
 | 
					    });
 | 
				
			||||||
#endif
 | 
					#endif
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
  
 | 
					  
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -110,12 +110,19 @@ int main (int argc, char** argv)
 | 
				
			|||||||
    PokeIndex<LorentzIndex>(mom, mommu, mu);
 | 
					    PokeIndex<LorentzIndex>(mom, mommu, mu);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    // fourth order exponential approx
 | 
					    // fourth order exponential approx
 | 
				
			||||||
    thread_loop( (auto i=mom.begin(); i<mom.end(); i++),{
 | 
					    auto mom_v    = mom.View();
 | 
				
			||||||
      Uprime[i](mu) = U[i](mu) + mom[i](mu)*U[i](mu)*dt + mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
 | 
					    auto U_v      = U.View();
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
 | 
					    auto Uprime_v = Uprime.View();
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
 | 
					
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
					    thread_loop( (auto i=mom_v.begin();i<mom_v.end();i++),{
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0);
 | 
					      Uprime_v[i](mu) =	  U_v[i](mu)
 | 
				
			||||||
 | 
						+ mom_v[i](mu)*U_v[i](mu)*dt 
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
				
			||||||
 | 
						;
 | 
				
			||||||
    });
 | 
					    });
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -120,15 +120,18 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
					    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    // fourth order exponential approx
 | 
					    // fourth order exponential approx
 | 
				
			||||||
    thread_loop( (auto i=mom.begin();i<mom.end();i++),{
 | 
					    auto mom_v    = mom.View();
 | 
				
			||||||
      Uprime[i](mu) =
 | 
					    auto U_v      = U.View();
 | 
				
			||||||
	  U[i](mu)
 | 
					    auto Uprime_v = Uprime.View();
 | 
				
			||||||
	+ mom[i](mu)*U[i](mu)*dt 
 | 
					
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
 | 
					    thread_loop( (auto i=mom_v.begin();i<mom_v.end();i++),{
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
 | 
					      Uprime_v[i](mu) =	  U_v[i](mu)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
 | 
						+ mom_v[i](mu)*U_v[i](mu)*dt 
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
				
			||||||
	;
 | 
						;
 | 
				
			||||||
    });
 | 
					    });
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -31,7 +31,6 @@ See the full license in the file "LICENSE" in the top level distribution directo
 | 
				
			|||||||
 | 
					
 | 
				
			||||||
using namespace std;
 | 
					using namespace std;
 | 
				
			||||||
using namespace Grid;
 | 
					using namespace Grid;
 | 
				
			||||||
 ;
 | 
					 | 
				
			||||||
 | 
					
 | 
				
			||||||
typedef GparityWilsonImplR FermionImplPolicy;
 | 
					typedef GparityWilsonImplR FermionImplPolicy;
 | 
				
			||||||
typedef GparityDomainWallEOFAFermionR FermionAction;
 | 
					typedef GparityDomainWallEOFAFermionR FermionAction;
 | 
				
			||||||
@@ -115,12 +114,19 @@ int main (int argc, char** argv)
 | 
				
			|||||||
    PokeIndex<LorentzIndex>(mom, mommu, mu);
 | 
					    PokeIndex<LorentzIndex>(mom, mommu, mu);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    // fourth order exponential approx
 | 
					    // fourth order exponential approx
 | 
				
			||||||
    thread_loop( (auto i=mom.begin(); i<mom.end(); i++),{
 | 
					    auto mom_v    = mom.View();
 | 
				
			||||||
      Uprime[i](mu) = U[i](mu) + mom[i](mu)*U[i](mu)*dt + mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
 | 
					    auto U_v      = U.View();
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
 | 
					    auto Uprime_v = Uprime.View();
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
 | 
					
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
					    thread_loop( (auto i=mom_v.begin();i<mom_v.end();i++),{
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0);
 | 
					      Uprime_v[i](mu) =	  U_v[i](mu)
 | 
				
			||||||
 | 
						+ mom_v[i](mu)*U_v[i](mu)*dt 
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
				
			||||||
 | 
						;
 | 
				
			||||||
    });
 | 
					    });
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -87,8 +87,11 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
					    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    // fourth order exponential approx
 | 
					    // fourth order exponential approx
 | 
				
			||||||
    thread_loop( (auto i=mom.begin();i<mom.end();i++),{ // exp(pmu dt) * Umu
 | 
					    auto mom_v   = mom.View();
 | 
				
			||||||
      Uprime[i](mu) = U[i](mu) + mom[i](mu)*U[i](mu)*dt ;
 | 
					    auto Uprime_v= Uprime.View();
 | 
				
			||||||
 | 
					    auto U_v     = U.View();
 | 
				
			||||||
 | 
					    thread_loop( (auto i=mom_v.begin();i<mom_v.end();i++),{ // exp(pmu dt) * Umu
 | 
				
			||||||
 | 
					      Uprime_v[i](mu) = U_v[i](mu) + mom_v[i](mu)*U_v[i](mu)*dt ;
 | 
				
			||||||
    });
 | 
					    });
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -105,15 +105,17 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
					    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    // fourth order exponential approx
 | 
					    // fourth order exponential approx
 | 
				
			||||||
    thread_loop( (auto i=mom.begin();i<mom.end();i++),{
 | 
					    auto U_v = U.View();
 | 
				
			||||||
      Uprime[i](mu) =
 | 
					    auto mom_v = mom.View();
 | 
				
			||||||
	  U[i](mu)
 | 
					    auto Uprime_v = Uprime.View();
 | 
				
			||||||
	+ mom[i](mu)*U[i](mu)*dt 
 | 
					    thread_loop( (auto i=mom_v.begin();i<mom_v.end();i++),{
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
 | 
					      Uprime_v[i](mu) = U_v[i](mu)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
 | 
						+ mom_v[i](mu)*U_v[i](mu)*dt 
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
				
			||||||
	;
 | 
						;
 | 
				
			||||||
    });
 | 
					    });
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -106,17 +106,21 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
					    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    // fourth order exponential approx
 | 
					    // fourth order exponential approx
 | 
				
			||||||
    thread_loop( (auto i=mom.begin();i<mom.end();i++),{
 | 
					    auto mom_v    = mom.View();
 | 
				
			||||||
      Uprime[i](mu) =
 | 
					    auto U_v      = U.View();
 | 
				
			||||||
	  U[i](mu)
 | 
					    auto Uprime_v = Uprime.View();
 | 
				
			||||||
	+ mom[i](mu)*U[i](mu)*dt 
 | 
					
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
 | 
					    thread_loop( (auto i=mom_v.begin();i<mom_v.end();i++),{
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
 | 
					      Uprime_v[i](mu) =	  U_v[i](mu)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
 | 
						+ mom_v[i](mu)*U_v[i](mu)*dt 
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
				
			||||||
	;
 | 
						;
 | 
				
			||||||
    });
 | 
					    });
 | 
				
			||||||
 | 
					
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
  
 | 
					  
 | 
				
			||||||
  Wil.ImportGauge(Uprime);
 | 
					  Wil.ImportGauge(Uprime);
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -101,15 +101,17 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
					    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    // fourth order exponential approx
 | 
					    // fourth order exponential approx
 | 
				
			||||||
    thread_loop( (auto i=mom.begin();i<mom.end();i++),{
 | 
					    auto U_v      = U.View();
 | 
				
			||||||
      Uprime[i](mu) =
 | 
					    auto mom_v    = mom.View();
 | 
				
			||||||
	  U[i](mu)
 | 
					    auto Uprime_v = Uprime.View();
 | 
				
			||||||
	+ mom[i](mu)*U[i](mu)*dt 
 | 
					    thread_loop( (auto i=mom_v.begin();i<mom_v.end();i++),{
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
 | 
					      Uprime_v[i](mu) =  U_v[i](mu)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
 | 
						+ mom_v[i](mu)*U_v[i](mu)*dt 
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
				
			||||||
	;
 | 
						;
 | 
				
			||||||
    });
 | 
					    });
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -112,12 +112,19 @@ int main (int argc, char** argv)
 | 
				
			|||||||
    PokeIndex<LorentzIndex>(mom, mommu, mu);
 | 
					    PokeIndex<LorentzIndex>(mom, mommu, mu);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    // fourth order exponential approx
 | 
					    // fourth order exponential approx
 | 
				
			||||||
    thread_loop( (auto i=mom.begin(); i<mom.end(); i++),{
 | 
					    auto mom_v    = mom.View();
 | 
				
			||||||
      Uprime[i](mu) = U[i](mu) + mom[i](mu)*U[i](mu)*dt + mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
 | 
					    auto U_v      = U.View();
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
 | 
					    auto Uprime_v = Uprime.View();
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
 | 
					
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
					    thread_loop( (auto i=mom_v.begin();i<mom_v.end();i++),{
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0);
 | 
					      Uprime_v[i](mu) =	  U_v[i](mu)
 | 
				
			||||||
 | 
						+ mom_v[i](mu)*U_v[i](mu)*dt 
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
				
			||||||
 | 
						;
 | 
				
			||||||
    });
 | 
					    });
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -115,14 +115,16 @@ int main (int argc, char** argv)
 | 
				
			|||||||
    SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
 | 
					    SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    PokeIndex<LorentzIndex>(mom, mommu, mu);
 | 
					    PokeIndex<LorentzIndex>(mom, mommu, mu);
 | 
				
			||||||
 | 
					    auto U_v      = U.View();
 | 
				
			||||||
 | 
					    auto mom_v    = mom.View();
 | 
				
			||||||
 | 
					    auto Uprime_v = Uprime.View();
 | 
				
			||||||
    // fourth order exponential approx
 | 
					    // fourth order exponential approx
 | 
				
			||||||
    thread_loop( (auto i=mom.begin(); i<mom.end(); i++),{
 | 
					    thread_loop( (auto i=mom_v.begin(); i<mom_v.end(); i++),{
 | 
				
			||||||
      Uprime[i](mu) = U[i](mu) + mom[i](mu)*U[i](mu)*dt + mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
 | 
					      Uprime_v[i](mu) = U_v[i](mu) + mom_v[i](mu)*U_v[i](mu)*dt + mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0)
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
 | 
					                        + mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0)
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
 | 
					                        + mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0)
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
					                        + mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
				
			||||||
                        + mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0);
 | 
					                        + mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0);
 | 
				
			||||||
    });
 | 
					    });
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -101,15 +101,17 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
					    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    // fourth order exponential approx
 | 
					    // fourth order exponential approx
 | 
				
			||||||
    thread_loop( (auto i=mom.begin();i<mom.end();i++),{
 | 
					    auto U_v = U.View();
 | 
				
			||||||
      Uprime[i](mu) =
 | 
					    auto mom_v = mom.View();
 | 
				
			||||||
	  U[i](mu)
 | 
					    auto Uprime_v = Uprime.View();
 | 
				
			||||||
	+ mom[i](mu)*U[i](mu)*dt 
 | 
					    thread_loop( (auto i=mom_v.begin();i<mom_v.end();i++),{
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
 | 
					      Uprime_v[i](mu) =  U_v[i](mu)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
 | 
						+ mom_v[i](mu)*U_v[i](mu)*dt 
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
				
			||||||
	;
 | 
						;
 | 
				
			||||||
    });
 | 
					    });
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -86,8 +86,11 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
					    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    // fourth order exponential approx
 | 
					    // fourth order exponential approx
 | 
				
			||||||
    thread_loop( (auto i=mom.begin();i<mom.end();i++),{ // exp(pmu dt) * Umu
 | 
					    auto Uprime_v = Uprime.View();
 | 
				
			||||||
      Uprime[i](mu) = U[i](mu) + mom[i](mu)*U[i](mu)*dt ;
 | 
					    auto U_v      = U.View();
 | 
				
			||||||
 | 
					    auto mom_v    = mom.View();
 | 
				
			||||||
 | 
					    thread_loop( (auto i=mom_v.begin();i<mom_v.end();i++),{ // exp(pmu dt) * Umu
 | 
				
			||||||
 | 
					      Uprime_v[i](mu) = U_v[i](mu) + mom_v[i](mu)*U_v[i](mu)*dt ;
 | 
				
			||||||
    });
 | 
					    });
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -29,9 +29,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 | 
				
			|||||||
 | 
					
 | 
				
			||||||
using namespace std;
 | 
					using namespace std;
 | 
				
			||||||
using namespace Grid;
 | 
					using namespace Grid;
 | 
				
			||||||
 ;
 | 
					 | 
				
			||||||
 | 
					 | 
				
			||||||
 
 | 
					 | 
				
			||||||
 | 
					
 | 
				
			||||||
int main (int argc, char ** argv)
 | 
					int main (int argc, char ** argv)
 | 
				
			||||||
{
 | 
					{
 | 
				
			||||||
@@ -104,14 +101,17 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
					    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    // fourth order exponential approx
 | 
					    // fourth order exponential approx
 | 
				
			||||||
    thread_loop( (auto i=mom.begin();i<mom.end();i++) ,{
 | 
					    auto U_v = U.View();
 | 
				
			||||||
      Uprime[i](mu)  =	  U[i](mu);
 | 
					    auto mom_v = mom.View();
 | 
				
			||||||
      Uprime[i](mu) += mom[i](mu)*U[i](mu)*dt ;
 | 
					    auto Uprime_v = Uprime.View();
 | 
				
			||||||
      Uprime[i](mu) += mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0);
 | 
					    thread_loop( (auto i=mom_v.begin();i<mom_v.end();i++) ,{
 | 
				
			||||||
      Uprime[i](mu) += mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0);
 | 
					      Uprime_v[i](mu)  =	  U_v[i](mu);
 | 
				
			||||||
      Uprime[i](mu) += mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0);
 | 
					      Uprime_v[i](mu) += mom_v[i](mu)*U_v[i](mu)*dt ;
 | 
				
			||||||
      Uprime[i](mu) += mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0);
 | 
					      Uprime_v[i](mu) += mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0);
 | 
				
			||||||
      Uprime[i](mu) += mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0);
 | 
					      Uprime_v[i](mu) += mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0);
 | 
				
			||||||
 | 
					      Uprime_v[i](mu) += mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0);
 | 
				
			||||||
 | 
					      Uprime_v[i](mu) += mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0);
 | 
				
			||||||
 | 
					      Uprime_v[i](mu) += mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0);
 | 
				
			||||||
    });
 | 
					    });
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 | 
					
 | 
				
			||||||
@@ -125,7 +125,6 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
  // Use derivative to estimate dS
 | 
					  // Use derivative to estimate dS
 | 
				
			||||||
  //////////////////////////////////////////////
 | 
					  //////////////////////////////////////////////
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 | 
					 | 
				
			||||||
  for(int mu=0;mu<Nd;mu++){
 | 
					  for(int mu=0;mu<Nd;mu++){
 | 
				
			||||||
    std::cout << "" <<std::endl;
 | 
					    std::cout << "" <<std::endl;
 | 
				
			||||||
    mommu   = PeekIndex<LorentzIndex>(mom,mu);
 | 
					    mommu   = PeekIndex<LorentzIndex>(mom,mu);
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -115,15 +115,18 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
					    PokeIndex<LorentzIndex>(mom,mommu,mu);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    // fourth order exponential approx
 | 
					    // fourth order exponential approx
 | 
				
			||||||
    thread_loop( (auto i=mom.begin();i<mom.end();i++),{
 | 
					    auto mom_v    = mom.View();
 | 
				
			||||||
      Uprime[i](mu) =
 | 
					    auto U_v      = U.View();
 | 
				
			||||||
	  U[i](mu)
 | 
					    auto Uprime_v = Uprime.View();
 | 
				
			||||||
	+ mom[i](mu)*U[i](mu)*dt 
 | 
					
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
 | 
					    thread_loop( (auto i=mom_v.begin();i<mom_v.end();i++),{
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
 | 
					      Uprime_v[i](mu) =	  U_v[i](mu)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
 | 
						+ mom_v[i](mu)*U_v[i](mu)*dt 
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0)
 | 
				
			||||||
	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0)
 | 
				
			||||||
 | 
						+ mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
 | 
				
			||||||
	;
 | 
						;
 | 
				
			||||||
    });
 | 
					    });
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -30,9 +30,9 @@ directory
 | 
				
			|||||||
/*  END LEGAL */
 | 
					/*  END LEGAL */
 | 
				
			||||||
#include <Grid/Grid.h>
 | 
					#include <Grid/Grid.h>
 | 
				
			||||||
 | 
					
 | 
				
			||||||
int main(int argc, char **argv) {
 | 
					int main(int argc, char **argv) 
 | 
				
			||||||
 | 
					{
 | 
				
			||||||
  using namespace Grid;
 | 
					  using namespace Grid;
 | 
				
			||||||
   ;
 | 
					 | 
				
			||||||
 | 
					
 | 
				
			||||||
  Grid_init(&argc, &argv);
 | 
					  Grid_init(&argc, &argv);
 | 
				
			||||||
  GridLogLayout();
 | 
					  GridLogLayout();
 | 
				
			||||||
 
 | 
				
			|||||||
		Reference in New Issue
	
	Block a user