portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

Big commit fixing nocompiles in defective C++11 compilers (gcc, icpc). stared getting to

Browse Source 
near the bleeding edge I guess
This commit is contained in:
Peter Boyle 2015-06-30 15:17:27 +01:00
parent 98c817df1b
commit 03ca506a3d
16 changed files with 349 additions and 101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
lib/Simd.h
View File

@ -17,6 +17,7 @@ namespace Grid {
typedef float RealF;
typedef double RealD;
#define GRID_DEFAULT_PRECISION_DOUBLE
#ifdef GRID_DEFAULT_PRECISION_DOUBLE
typedef RealD Real;
#else

7
lib/algorithms/CoarsenedMatrix.h
View File

@ -121,7 +121,7 @@ namespace Grid {
RealD scale;
ConjugateGradient<FineField> CG(1.0e-3,10000);
ConjugateGradient<FineField> CG(1.0e-4,10000);
FineField noise(FineGrid);
FineField Mn(FineGrid);
@ -131,7 +131,8 @@ namespace Grid {
scale = std::pow(norm2(noise),-0.5);
noise=noise*scale;
hermop.Op(noise,Mn); std::cout << "Noise "<<b<<" nMdagMn "<<norm2(Mn)<< " "<< norm2(noise)<<std::endl;
hermop.Op(noise,Mn); std::cout << "noise ["<<b<<"] <n|MdagM|n> "<<norm2(Mn)<<std::endl;
for(int i=0;i<2;i++){
CG(hermop,noise,subspace[b]);
@ -142,7 +143,7 @@ namespace Grid {
}
hermop.Op(noise,Mn); std::cout << "filtered "<<b<<" <s|MdagM|s> "<<norm2(Mn)<< " "<< norm2(noise)<<std::endl;
hermop.Op(noise,Mn); std::cout << "filtered["<<b<<"] <f|MdagM|f> "<<norm2(Mn)<<std::endl;
subspace[b] = noise;
}

26
lib/lattice/Lattice_peekpoke.h
View File

@ -10,18 +10,8 @@ namespace Grid {
////////////////////////////////////////////////////////////////////////////////////////////////////
// Peek internal indices of a Lattice object
////////////////////////////////////////////////////////////////////////////////////////////////////
template<int Index,class vobj>
auto peekIndex(const Lattice<vobj> &lhs) -> Lattice<decltype(peekIndex<Index>(lhs._odata[0]))>
{
Lattice<decltype(peekIndex<Index>(lhs._odata[0]))> ret(lhs._grid);
PARALLEL_FOR_LOOP
for(int ss=0;ss<lhs._grid->oSites();ss++){
ret._odata[ss] = peekIndex<Index>(lhs._odata[ss]);
}
return ret;
};
template<int Index,class vobj>
auto peekIndex(const Lattice<vobj> &lhs,int i) -> Lattice<decltype(peekIndex<Index>(lhs._odata[0],i))>
auto PeekIndex(const Lattice<vobj> &lhs,int i) -> Lattice<decltype(peekIndex<Index>(lhs._odata[0],i))>
{
Lattice<decltype(peekIndex<Index>(lhs._odata[0],i))> ret(lhs._grid);
PARALLEL_FOR_LOOP
@ -31,7 +21,7 @@ PARALLEL_FOR_LOOP
return ret;
};
template<int Index,class vobj>
auto peekIndex(const Lattice<vobj> &lhs,int i,int j) -> Lattice<decltype(peekIndex<Index>(lhs._odata[0],i,j))>
auto PeekIndex(const Lattice<vobj> &lhs,int i,int j) -> Lattice<decltype(peekIndex<Index>(lhs._odata[0],i,j))>
{
Lattice<decltype(peekIndex<Index>(lhs._odata[0],i,j))> ret(lhs._grid);
PARALLEL_FOR_LOOP
@ -44,16 +34,8 @@ PARALLEL_FOR_LOOP
////////////////////////////////////////////////////////////////////////////////////////////////////
// Poke internal indices of a Lattice object
////////////////////////////////////////////////////////////////////////////////////////////////////
template<int Index,class vobj>
void pokeIndex(Lattice<vobj> &lhs,const Lattice<decltype(peekIndex<Index>(lhs._odata[0]))> & rhs)
{
PARALLEL_FOR_LOOP
for(int ss=0;ss<lhs._grid->oSites();ss++){
pokeIndex<Index>(lhs._odata[ss],rhs._odata[ss]);
}
}
template<int Index,class vobj>
void pokeIndex(Lattice<vobj> &lhs,const Lattice<decltype(peekIndex<Index>(lhs._odata[0],0))> & rhs,int i)
void PokeIndex(Lattice<vobj> &lhs,const Lattice<decltype(peekIndex<Index>(lhs._odata[0],0))> & rhs,int i)
{
PARALLEL_FOR_LOOP
for(int ss=0;ss<lhs._grid->oSites();ss++){
@ -61,7 +43,7 @@ PARALLEL_FOR_LOOP
}
}
template<int Index,class vobj>
void pokeIndex(Lattice<vobj> &lhs,const Lattice<decltype(peekIndex<Index>(lhs._odata[0],0,0))> & rhs,int i,int j)
void PokeIndex(Lattice<vobj> &lhs,const Lattice<decltype(peekIndex<Index>(lhs._odata[0],0,0))> & rhs,int i,int j)
{
PARALLEL_FOR_LOOP
for(int ss=0;ss<lhs._grid->oSites();ss++){

2
lib/lattice/Lattice_trace.h
View File

@ -26,7 +26,7 @@ PARALLEL_FOR_LOOP
// Trace Index level dependent operation
////////////////////////////////////////////////////////////////////////////////////////////////////
template<int Index,class vobj>
inline auto traceIndex(const Lattice<vobj> &lhs) -> Lattice<decltype(traceIndex<Index>(lhs._odata[0]))>
inline auto TraceIndex(const Lattice<vobj> &lhs) -> Lattice<decltype(traceIndex<Index>(lhs._odata[0]))>
{
Lattice<decltype(traceIndex<Index>(lhs._odata[0]))> ret(lhs._grid);
PARALLEL_FOR_LOOP

2
lib/lattice/Lattice_transpose.h
View File

@ -24,7 +24,7 @@ PARALLEL_FOR_LOOP
// Index level dependent transpose
////////////////////////////////////////////////////////////////////////////////////////////////////
template<int Index,class vobj>
inline auto transposeIndex(const Lattice<vobj> &lhs) -> Lattice<decltype(transposeIndex<Index>(lhs._odata[0]))>
inline auto TransposeIndex(const Lattice<vobj> &lhs) -> Lattice<decltype(transposeIndex<Index>(lhs._odata[0]))>
{
Lattice<decltype(transposeIndex<Index>(lhs._odata[0]))> ret(lhs._grid);
PARALLEL_FOR_LOOP

4
lib/qcd/action/fermion/WilsonFermion.cc
View File

@ -36,9 +36,9 @@ void WilsonFermion::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGauge
LatticeColourMatrix U(GaugeGrid());
for(int mu=0;mu<Nd;mu++){
U = PeekIndex<LorentzIndex>(Umu,mu);
pokeIndex<LorentzIndex>(Uds,U,mu);
PokeIndex<LorentzIndex>(Uds,U,mu);
U = adj(Cshift(U,mu,-1));
pokeIndex<LorentzIndex>(Uds,U,mu+4);
PokeIndex<LorentzIndex>(Uds,U,mu+4);
}
}

4
lib/qcd/action/fermion/WilsonFermion5D.cc
View File

@ -78,9 +78,9 @@ void WilsonFermion5D::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGau
LatticeColourMatrix U(GaugeGrid());
for(int mu=0;mu<Nd;mu++){
U = PeekIndex<LorentzIndex>(Umu,mu);
pokeIndex<LorentzIndex>(Uds,U,mu);
PokeIndex<LorentzIndex>(Uds,U,mu);
U = adj(Cshift(U,mu,-1));
pokeIndex<LorentzIndex>(Uds,U,mu+4);
PokeIndex<LorentzIndex>(Uds,U,mu+4);
}
}
void WilsonFermion5D::DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir5,int disp)

6
lib/qcd/utils/SUn.h
View File

@ -547,21 +547,21 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
LatticeMatrix Umu(out._grid);
for(int mu=0;mu<Nd;mu++){
LieRandomize(pRNG,Umu,1.0);
pokeIndex<LorentzIndex>(out,Umu,mu);
PokeIndex<LorentzIndex>(out,Umu,mu);
}
}
static void TepidConfiguration(GridParallelRNG &pRNG,LatticeGaugeField &out){
LatticeMatrix Umu(out._grid);
for(int mu=0;mu<Nd;mu++){
LieRandomize(pRNG,Umu,0.01);
pokeLorentz(out,Umu,mu);
PokeLorentz(out,Umu,mu);
}
}
static void ColdConfiguration(GridParallelRNG &pRNG,LatticeGaugeField &out){
LatticeMatrix Umu(out._grid);
Umu=1.0;
for(int mu=0;mu<Nd;mu++){
pokeLorentz(out,Umu,mu);
PokeLorentz(out,Umu,mu);
}
}

18
lib/tensors/Tensor_class.h
View File

@ -91,18 +91,12 @@ public:
return _internal;
}
// Type casts meta programmed
template<class U=vtype,class V=scalar_type,IfComplex<V> = 0,IfNotSimd<U> = 0>
operator ComplexF () const { return(TensorRemove(_internal)); };
template<class U=vtype,class V=scalar_type,IfComplex<V> = 0,IfNotSimd<U> = 0>
operator ComplexD () const { return(TensorRemove(_internal)); };
template<class U=vtype,class V=scalar_type,IfComplex<V> = 0,IfNotSimd<U> = 0>
operator RealD () const { return(real(TensorRemove(_internal))); }
template<class U=vtype,class V=scalar_type,IfReal<V> = 0,IfNotSimd<U> = 0>
operator RealD () const { return TensorRemove(_internal); }
template<class U=vtype,class V=scalar_type,IfInteger<V> = 0,IfNotSimd<U> = 0>
operator Integer () const { return Integer(TensorRemove(_internal)); }
// Type casts meta programmed, must be pure scalar to match TensorRemove
template<class U=vtype,class V=scalar_type,IfComplex<V> = 0,IfNotSimd<U> = 0> operator ComplexF () const { return(TensorRemove(_internal)); };
template<class U=vtype,class V=scalar_type,IfComplex<V> = 0,IfNotSimd<U> = 0> operator ComplexD () const { return(TensorRemove(_internal)); };
// template<class U=vtype,class V=scalar_type,IfComplex<V> = 0,IfNotSimd<U> = 0> operator RealD () const { return(real(TensorRemove(_internal))); }
template<class U=vtype,class V=scalar_type,IfReal<V> = 0,IfNotSimd<U> = 0> operator RealD () const { return TensorRemove(_internal); }
template<class U=vtype,class V=scalar_type,IfInteger<V> = 0,IfNotSimd<U> = 0> operator Integer () const { return Integer(TensorRemove(_internal)); }
// convert from a something to a scalar via constructor of something arg
template<class T,typename std::enable_if<!isGridTensor<T>::value, T>::type* = nullptr > strong_inline iScalar<vtype> operator = (T arg)

4
lib/tensors/Tensor_peek.h
View File

@ -11,7 +11,7 @@ namespace Grid {
//template<int Level> inline ComplexD peekIndex(const ComplexD arg) { return arg;}
//template<int Level> inline RealF peekIndex(const RealF arg) { return arg;}
//template<int Level> inline RealD peekIndex(const RealD arg) { return arg;}
#if 0
// Scalar peek, no indices
template<int Level,class vtype,typename std::enable_if< iScalar<vtype>::TensorLevel == Level >::type * =nullptr> inline
auto peekIndex(const iScalar<vtype> &arg) -> iScalar<vtype>
@ -88,6 +88,7 @@ template<int Level,class vtype,int N, typename std::enable_if< iScalar<vtype>::T
}
return ret;
}
// matrix
template<int Level,class vtype,int N, typename std::enable_if< iScalar<vtype>::TensorLevel != Level >::type * =nullptr> inline
auto peekIndex(const iMatrix<vtype,N> &arg) -> iMatrix<decltype(peekIndex<Level>(arg._internal[0][0])),N>
@ -119,6 +120,7 @@ template<int Level,class vtype,int N, typename std::enable_if< iScalar<vtype>::T
}}
return ret;
}
#endif
}

7
lib/tensors/Tensor_poke.h
View File

@ -5,7 +5,7 @@ namespace Grid {
//////////////////////////////////////////////////////////////////////////////
// Poke a specific index;
//////////////////////////////////////////////////////////////////////////////
#if 0
// Scalar poke
template<int Level,class vtype,typename std::enable_if< iScalar<vtype>::TensorLevel == Level >::type * =nullptr> inline
void pokeIndex(iScalar<vtype> &ret, const iScalar<vtype> &arg)
@ -18,7 +18,7 @@ template<int Level,class vtype,int N,typename std::enable_if< iScalar<vtype>::Te
{
ret._internal[i] = arg._internal;
}
// Vector poke, two indices
//Matrix poke, two indices
template<int Level,class vtype,int N,typename std::enable_if< iScalar<vtype>::TensorLevel == Level >::type * =nullptr> inline
void pokeIndex(iMatrix<vtype,N> &ret, const iScalar<vtype> &arg,int i,int j)
{
@ -31,7 +31,6 @@ template<int Level,class vtype,int N,typename std::enable_if< iScalar<vtype>::Te
// scalar
template<int Level,class vtype,typename std::enable_if< iScalar<vtype>::TensorLevel != Level >::type * =nullptr> inline
void pokeIndex(iScalar<vtype> &ret, const iScalar<decltype(peekIndex<Level>(ret._internal))> &arg)
{
pokeIndex<Level>(ret._internal,arg._internal);
}
@ -95,7 +94,7 @@ template<int Level,class vtype,int N,typename std::enable_if< iScalar<vtype>::Te
pokeIndex<Level>(ret._internal[ii][jj],arg._internal[ii][jj],i,j);
}}
}
#endif
}
#endif

361
lib/tensors/Tensor_trace.h
View File

@ -32,58 +32,327 @@ inline auto trace(const iScalar<vtype> &arg) -> iScalar<decltype(trace(arg._inte
ret._internal=trace(arg._internal);
return ret;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////
// Trace Specific indices.
////////////////////////////////////////////////////////////////////////////////////////////////////////
template<int Level,class vtype,typename std::enable_if< iScalar<vtype>::TensorLevel != Level >::type * =nullptr> inline auto
traceIndex (const iScalar<vtype> &arg) -> iScalar<decltype(traceIndex<Level>(arg._internal))>
{
iScalar<decltype(traceIndex<Level>(arg._internal))> ret;
ret._internal=traceIndex<Level>(arg._internal);
return ret;
}
template<int Level,class vtype,typename std::enable_if< iScalar<vtype>::TensorLevel == Level >::type * =nullptr> inline auto
traceIndex (const iScalar<vtype> &arg) -> iScalar<vtype>
{
return arg;
}
// If we hit the right index, return scalar and trace it with no further recursion
template<int Level,class vtype,int N,typename std::enable_if< iScalar<vtype>::TensorLevel == Level >::type * =nullptr> inline
auto traceIndex(const iMatrix<vtype,N> &arg) -> iScalar<vtype>
{
iScalar<vtype> ret;
zeroit(ret._internal);
for(int i=0;i<N;i++){
ret._internal = ret._internal + arg._internal[i][i];
}
return ret;
}
// not this level, so recurse
template<int Level,class vtype,int N,typename std::enable_if< iScalar<vtype>::TensorLevel != Level >::type * =nullptr> inline
auto traceIndex(const iMatrix<vtype,N> &arg) -> iMatrix<decltype(traceIndex<Level>(arg._internal[0][0])),N>
{
iMatrix<decltype(traceIndex<Level>(arg._internal[0][0])),N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j] = traceIndex<Level>(arg._internal[i][j]);
}}
return ret;
}
/////////////////////////////////////
// Alternate implementation .. count
// indices from left. Annoying but
// I'm working around compiler bugs in gcc and earlier icpc
/////////////////////////////////////
// Allow to recurse if vector, but never terminate on a vector
// trace of a different index can distribute across the vector index in a replicated way
// but we do not trace a vector index.
template<int Level,class vtype,int N,typename std::enable_if< iVector<vtype, N>::TensorLevel != Level >::type * =nullptr> inline
auto traceIndex(const iVector<vtype,N> &arg) -> iVector<decltype(traceIndex<Level>(arg._internal[0])),N>
{
iVector<decltype(traceIndex<Level>(arg._internal[0])),N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = traceIndex<Level>(arg._internal[i]);
// ltrace of a different index can distribute across the vector index in a replicated way
// but we do not ltrace a vector index.
template<int Level>
class TensorIndexRecursion {
public:
////////////////////////////////////////////
// Recursion for tracing a specific index
////////////////////////////////////////////
template<class vtype>
static auto traceIndex(const iScalar<vtype> arg) -> iScalar<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal))>
{
iScalar<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal))> ret;
ret._internal = TensorIndexRecursion<Level-1>::traceIndex(arg._internal);
return ret;
}
template<class vtype,int N>
static auto traceIndex(const iVector<vtype,N> arg) -> iVector<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal[0])),N>
{
iVector<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal[0])),N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = TensorIndexRecursion<Level-1>::traceIndex(arg._internal[i]);
}
return ret;
}
template<class vtype,int N>
static auto traceIndex(const iMatrix<vtype,N> arg) -> iMatrix<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal[0][0])),N>
{
iMatrix<decltype(TensorIndexRecursion<Level-1>::traceIndex(arg._internal[0][0])),N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j] = TensorIndexRecursion<Level-1>::traceIndex(arg._internal[i][j]);
}}
return ret;
}
////////////////////////////////////////////
// Recursion for peeking a specific index
////////////////////////////////////////////
template<class vtype>
static auto peekIndex(const iScalar<vtype> arg,int i) -> iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal,0))>
{
iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal,0))> ret;
ret._internal = TensorIndexRecursion<Level-1>::peekIndex(arg._internal,i);
return ret;
}
template<class vtype>
static auto peekIndex(const iScalar<vtype> arg,int i,int j) -> iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal,0,0))>
{
iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal,0,0))> ret;
ret._internal = TensorIndexRecursion<Level-1>::peekIndex(arg._internal,i,j);
return ret;
}
template<class vtype,int N>
static auto peekIndex(const iVector<vtype,N> arg,int ii) -> iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0],0)),N>
{
iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0],0)),N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = TensorIndexRecursion<Level-1>::peekIndex(arg._internal[i],ii);
}
return ret;
}
template<class vtype,int N>
static auto peekIndex(const iVector<vtype,N> arg,int ii,int jj) -> iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0],0,0)),N>
{
iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0],0,0)),N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = TensorIndexRecursion<Level-1>::peekIndex(arg._internal[i],ii,jj);
}
return ret;
}
template<class vtype,int N>
static auto peekIndex(const iMatrix<vtype,N> arg,int ii) -> iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0][0],0)),N>
{
iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0][0],0)),N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j] = TensorIndexRecursion<Level-1>::peekIndex(arg._internal[i][j],ii);
}}
return ret;
}
template<class vtype,int N>
static auto peekIndex(const iMatrix<vtype,N> arg,int ii,int jj) -> iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0][0],0,0)),N>
{
iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(arg._internal[0][0],0,0)),N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j] = TensorIndexRecursion<Level-1>::peekIndex(arg._internal[i][j],ii,jj);
}}
return ret;
}
////////////////////////////////////////////
// Recursion for poking a specific index
////////////////////////////////////////////
template<class vtype> inline static
void pokeIndex(iScalar<vtype> &ret, const iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0))> &arg, int i)
{
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal,arg._internal,i);
}
template<class vtype> inline static
void pokeIndex(iScalar<vtype> &ret, const iScalar<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0,0))> &arg, int i,int j)
{
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal,arg._internal,i,j);
}
template<class vtype,int N> inline static
void pokeIndex(iVector<vtype,N> &ret, const iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0)),N> &arg, int i)
{
for(int ii=0;ii<N;ii++){
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal[ii],arg._internal[ii],i);
}
}
template<class vtype,int N> inline static
void pokeIndex(iVector<vtype,N> &ret, const iVector<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0)),N> &arg, int i,int j)
{
for(int ii=0;ii<N;ii++){
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal[ii],arg._internal[ii],i,j);
}
}
template<class vtype,int N> inline static
void pokeIndex(iMatrix<vtype,N> &ret, const iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0)),N> &arg, int i)
{
for(int ii=0;ii<N;ii++){
for(int jj=0;jj<N;jj++){
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal[ii][jj],arg._internal[ii][jj],i);
}}
}
template<class vtype,int N> inline static
void pokeIndex(iMatrix<vtype,N> &ret, const iMatrix<decltype(TensorIndexRecursion<Level-1>::peekIndex(ret._internal,0)),N> &arg, int i,int j)
{
for(int ii=0;ii<N;ii++){
for(int jj=0;jj<N;jj++){
TensorIndexRecursion<Level-1>::pokeIndex(ret._internal[ii][jj],arg._internal[ii][jj],i,j);
}}
}
////////////////////////////////////////////
// Recursion for transposing a specific index
////////////////////////////////////////////
template<class vtype>
static auto transposeIndex(const iScalar<vtype> arg) -> iScalar<vtype>
{
iScalar<vtype> ret;
ret._internal = TensorIndexRecursion<Level-1>::transposeIndex(arg._internal);
return ret;
}
template<class vtype,int N>
static auto transposeIndex(const iVector<vtype,N> arg) -> iVector<vtype,N>
{
iVector<vtype,N> ret;
for(int i=0;i<N;i++){
ret._internal[i] = TensorIndexRecursion<Level-1>::transposeIndex(arg._internal[i]);
}
return ret;
}
template<class vtype,int N>
static auto transposeIndex(const iMatrix<vtype,N> arg) -> iMatrix<vtype,N>
{
iMatrix<vtype,N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j] = TensorIndexRecursion<Level-1>::transposeIndex(arg._internal[i][j]);
}}
return ret;
}
};
////////////////////////////
// strip const & ref quali's
////////////////////////////
#define RemoveCRV(a) typename std::remove_const<typename std::remove_reference<decltype(a)>::type>::type
template<>
class TensorIndexRecursion<0> {
public:
/////////////////////////////////////////
// Ends recursion for trace (scalar/vector/matrix)
/////////////////////////////////////////
template<class vtype>
static auto traceIndex(const iScalar<vtype> arg) -> iScalar<RemoveCRV(arg._internal)>
{
iScalar<RemoveCRV(arg._internal)> ret;
ret._internal = arg._internal;
return ret;
}
template<class vtype,int N>
static auto traceIndex(const iVector<vtype,N> arg) -> iScalar<RemoveCRV(arg._internal[0])>
{
iScalar<RemoveCRV(arg._internal[0])> ret;
ret._internal=zero;
for(int i=0;i<N;i++){
ret._internal = ret._internal+ arg._internal[i];
}
return ret;
}
template<class vtype,int N>
static auto traceIndex(const iMatrix<vtype,N> arg) -> iScalar<RemoveCRV(arg._internal[0][0])>
{
iScalar<RemoveCRV(arg._internal[0][0])> ret;
ret=zero;
for(int i=0;i<N;i++){
ret._internal = ret._internal+arg._internal[i][i];
}
return ret;
}
/////////////////////////////////////////
// Ends recursion for transpose scalar/vector/matrix
/////////////////////////////////////////
template<class vtype>
static auto transposeIndex(const iScalar<vtype> arg) -> iScalar<vtype>
{
iScalar<vtype> ret;
ret._internal = arg._internal;
return ret;
}
template<class vtype,int N>
static auto transposeIndex(const iVector<vtype,N> arg) -> iVector<vtype,N>
{
iScalar<vtype> ret;
ret._internal=zero;
for(int i=0;i<N;i++){
ret._internal = ret._internal+ arg._internal[i];
}
return ret;
}
template<class vtype,int N>
static auto transposeIndex(const iMatrix<vtype,N> arg) -> iMatrix<vtype,N>
{
iScalar<vtype> ret;
ret=zero;
for(int i=0;i<N;i++){
ret._internal = ret._internal+arg._internal[i][i];
}
return ret;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// End recursion for peeking a specific index; single index on vector, double index on matrix
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template<class vtype,int N>
static auto peekIndex(const iVector<vtype,N> arg,int ii) -> iScalar<vtype>
{
iScalar<vtype> ret;
ret._internal = arg._internal[ii];
return ret;
}
template<class vtype,int N>
static auto peekIndex(const iMatrix<vtype,N> arg,int ii,int jj) -> iScalar<vtype>
{
iScalar<vtype> ret;
ret._internal = arg._internal[ii][jj];
return ret;
}
// Vector poke, one index
template<class vtype,int N> inline static
void pokeIndex(iVector<vtype,N> &ret, const iScalar<vtype> &arg,int i)
{
ret._internal[i] = arg._internal;
}
// Matrix poke two indices
template<class vtype,int N> inline static
void pokeIndex(iMatrix<vtype,N> &ret, const iScalar<vtype> &arg,int i,int j)
{
ret._internal[i][j] = arg._internal;
}
};
////////////////////////////////////////////////////////////////////////////////////////////////////////
// External wrappers
////////////////////////////////////////////////////////////////////////////////////////////////////////
template<int Level,class vtype> inline auto traceIndex (const vtype &arg) -> RemoveCRV(TensorIndexRecursion<Level>::traceIndex(arg))
{
RemoveCRV(TensorIndexRecursion<Level>::traceIndex(arg)) ret;
ret=TensorIndexRecursion<Level>::traceIndex(arg);
return ret;
}
template<int Level,class vtype> inline auto transposeIndex (const vtype &arg) -> RemoveCRV(TensorIndexRecursion<Level>::transposeIndex(arg))
{
RemoveCRV(TensorIndexRecursion<Level>::transposeIndex(arg)) ret;
ret=TensorIndexRecursion<Level>::transposeIndex(arg);
return ret;
}
template<int Level,class vtype> inline auto peekIndex (const vtype &arg,int i) -> RemoveCRV(TensorIndexRecursion<Level>::peekIndex(arg,0))
{
RemoveCRV(TensorIndexRecursion<Level>::peekIndex(arg,0)) ret;
ret=TensorIndexRecursion<Level>::peekIndex(arg,i);
return ret;
}
template<int Level,class vtype> inline auto peekIndex (const vtype &arg,int i,int j) -> RemoveCRV(TensorIndexRecursion<Level>::peekIndex(arg,0,0))
{
RemoveCRV(TensorIndexRecursion<Level>::peekIndex(arg,0,0)) ret;
ret=TensorIndexRecursion<Level>::peekIndex(arg,i,j);
return ret;
}
template<int Level,class vtype> inline
void pokeIndex (vtype &ret,const decltype(TensorIndexRecursion<Level>::peekIndex(ret,0)) &arg,int i)
{
TensorIndexRecursion<Level>::pokeIndex(ret,arg,i);
}
template<int Level,class vtype> inline
void pokeIndex (vtype &ret,const decltype(TensorIndexRecursion<Level>::peekIndex(ret,0,0)) &arg,int i,int j)
{
TensorIndexRecursion<Level>::pokeIndex(ret,arg,i,j);
}
#undef RemoveCRV
}
#endif

2
tests/Test_dwf_even_odd.cc
View File

@ -53,7 +53,7 @@ int main (int argc, char ** argv)
random(RNG4,U[nn]);
if ( nn>0 )
U[nn]=zero;
pokeIndex<LorentzIndex>(Umu,U[nn],nn);
PokeIndex<LorentzIndex>(Umu,U[nn],nn);
}
RealD mass=0.1;

2
tests/Test_main.cc
View File

@ -141,7 +141,7 @@ int main (int argc, char ** argv)
// rscalar=real(scalar);
// iscalar=imag(scalar);
// scalar =cmplx(rscalar,iscalar);
pokeIndex<2>(cVec,scalar,1);
PokeIndex<ColourIndex>(cVec,scalar,1);
scalar=transpose(scalar);

2
tests/Test_quenched_update.cc
View File

@ -72,7 +72,7 @@ int main (int argc, char ** argv)
}
pokeIndex<LorentzIndex>(Umu,link,mu);
PokeIndex<LorentzIndex>(Umu,link,mu);
//reunitarise link;
}

2
tests/Test_wilson_even_odd.cc
View File

@ -55,7 +55,7 @@ int main (int argc, char ** argv)
Umu=zero;
for(int nn=0;nn<Nd;nn++){
random(pRNG,U[nn]);
pokeIndex<LorentzIndex>(Umu,U[nn],nn);
PokeIndex<LorentzIndex>(Umu,U[nn],nn);
}
RealD mass=0.1;