From f41c7dffef6ddeb643679afe9c383125b987923f Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 30 Jun 2015 15:17:27 +0100 Subject: [PATCH] Big commit fixing nocompiles in defective C++11 compilers (gcc, icpc). stared getting to near the bleeding edge I guess --- lib/Simd.h | 1 + lib/algorithms/CoarsenedMatrix.h | 7 +- lib/lattice/Lattice_peekpoke.h | 26 +- lib/lattice/Lattice_trace.h | 2 +- lib/lattice/Lattice_transpose.h | 2 +- lib/qcd/action/fermion/WilsonFermion.cc | 4 +- lib/qcd/action/fermion/WilsonFermion5D.cc | 4 +- lib/qcd/utils/SUn.h | 6 +- lib/tensors/Tensor_class.h | 18 +- lib/tensors/Tensor_peek.h | 4 +- lib/tensors/Tensor_poke.h | 7 +- lib/tensors/Tensor_trace.h | 361 +++++++++++++++++++--- tests/Test_dwf_even_odd.cc | 2 +- tests/Test_main.cc | 2 +- tests/Test_quenched_update.cc | 2 +- tests/Test_wilson_even_odd.cc | 2 +- 16 files changed, 349 insertions(+), 101 deletions(-) diff --git a/lib/Simd.h b/lib/Simd.h index 7da0d7ad..474e5dd3 100644 --- a/lib/Simd.h +++ b/lib/Simd.h @@ -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 diff --git a/lib/algorithms/CoarsenedMatrix.h b/lib/algorithms/CoarsenedMatrix.h index 9ea77b84..e257294c 100644 --- a/lib/algorithms/CoarsenedMatrix.h +++ b/lib/algorithms/CoarsenedMatrix.h @@ -121,7 +121,7 @@ namespace Grid { RealD scale; - ConjugateGradient CG(1.0e-3,10000); + ConjugateGradient 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 "< "< "< "< - auto peekIndex(const Lattice &lhs) -> Lattice(lhs._odata[0]))> - { - Lattice(lhs._odata[0]))> ret(lhs._grid); -PARALLEL_FOR_LOOP - for(int ss=0;ssoSites();ss++){ - ret._odata[ss] = peekIndex(lhs._odata[ss]); - } - return ret; - }; template - auto peekIndex(const Lattice &lhs,int i) -> Lattice(lhs._odata[0],i))> + auto PeekIndex(const Lattice &lhs,int i) -> Lattice(lhs._odata[0],i))> { Lattice(lhs._odata[0],i))> ret(lhs._grid); PARALLEL_FOR_LOOP @@ -31,7 +21,7 @@ PARALLEL_FOR_LOOP return ret; }; template - auto peekIndex(const Lattice &lhs,int i,int j) -> Lattice(lhs._odata[0],i,j))> + auto PeekIndex(const Lattice &lhs,int i,int j) -> Lattice(lhs._odata[0],i,j))> { Lattice(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 - void pokeIndex(Lattice &lhs,const Lattice(lhs._odata[0]))> & rhs) - { -PARALLEL_FOR_LOOP - for(int ss=0;ssoSites();ss++){ - pokeIndex(lhs._odata[ss],rhs._odata[ss]); - } - } template - void pokeIndex(Lattice &lhs,const Lattice(lhs._odata[0],0))> & rhs,int i) + void PokeIndex(Lattice &lhs,const Lattice(lhs._odata[0],0))> & rhs,int i) { PARALLEL_FOR_LOOP for(int ss=0;ssoSites();ss++){ @@ -61,7 +43,7 @@ PARALLEL_FOR_LOOP } } template - void pokeIndex(Lattice &lhs,const Lattice(lhs._odata[0],0,0))> & rhs,int i,int j) + void PokeIndex(Lattice &lhs,const Lattice(lhs._odata[0],0,0))> & rhs,int i,int j) { PARALLEL_FOR_LOOP for(int ss=0;ssoSites();ss++){ diff --git a/lib/lattice/Lattice_trace.h b/lib/lattice/Lattice_trace.h index 75cc5b87..1b35b412 100644 --- a/lib/lattice/Lattice_trace.h +++ b/lib/lattice/Lattice_trace.h @@ -26,7 +26,7 @@ PARALLEL_FOR_LOOP // Trace Index level dependent operation //////////////////////////////////////////////////////////////////////////////////////////////////// template - inline auto traceIndex(const Lattice &lhs) -> Lattice(lhs._odata[0]))> + inline auto TraceIndex(const Lattice &lhs) -> Lattice(lhs._odata[0]))> { Lattice(lhs._odata[0]))> ret(lhs._grid); PARALLEL_FOR_LOOP diff --git a/lib/lattice/Lattice_transpose.h b/lib/lattice/Lattice_transpose.h index 73ef1747..a23ce48c 100644 --- a/lib/lattice/Lattice_transpose.h +++ b/lib/lattice/Lattice_transpose.h @@ -24,7 +24,7 @@ PARALLEL_FOR_LOOP // Index level dependent transpose //////////////////////////////////////////////////////////////////////////////////////////////////// template - inline auto transposeIndex(const Lattice &lhs) -> Lattice(lhs._odata[0]))> + inline auto TransposeIndex(const Lattice &lhs) -> Lattice(lhs._odata[0]))> { Lattice(lhs._odata[0]))> ret(lhs._grid); PARALLEL_FOR_LOOP diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index 067d5113..f5ff7cb0 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -36,9 +36,9 @@ void WilsonFermion::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGauge LatticeColourMatrix U(GaugeGrid()); for(int mu=0;mu(Umu,mu); - pokeIndex(Uds,U,mu); + PokeIndex(Uds,U,mu); U = adj(Cshift(U,mu,-1)); - pokeIndex(Uds,U,mu+4); + PokeIndex(Uds,U,mu+4); } } diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 6d6097cf..316b0676 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -78,9 +78,9 @@ void WilsonFermion5D::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGau LatticeColourMatrix U(GaugeGrid()); for(int mu=0;mu(Umu,mu); - pokeIndex(Uds,U,mu); + PokeIndex(Uds,U,mu); U = adj(Cshift(U,mu,-1)); - pokeIndex(Uds,U,mu+4); + PokeIndex(Uds,U,mu+4); } } void WilsonFermion5D::DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir5,int disp) diff --git a/lib/qcd/utils/SUn.h b/lib/qcd/utils/SUn.h index 4088a23b..69de3f9f 100644 --- a/lib/qcd/utils/SUn.h +++ b/lib/qcd/utils/SUn.h @@ -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(out,Umu,mu); + PokeIndex(out,Umu,mu); } } static void TepidConfiguration(GridParallelRNG &pRNG,LatticeGaugeField &out){ LatticeMatrix Umu(out._grid); for(int mu=0;mu = 0,IfNotSimd = 0> - operator ComplexF () const { return(TensorRemove(_internal)); }; - template = 0,IfNotSimd = 0> - operator ComplexD () const { return(TensorRemove(_internal)); }; - template = 0,IfNotSimd = 0> - operator RealD () const { return(real(TensorRemove(_internal))); } - template = 0,IfNotSimd = 0> - operator RealD () const { return TensorRemove(_internal); } - template = 0,IfNotSimd = 0> - operator Integer () const { return Integer(TensorRemove(_internal)); } - + // Type casts meta programmed, must be pure scalar to match TensorRemove + template = 0,IfNotSimd = 0> operator ComplexF () const { return(TensorRemove(_internal)); }; + template = 0,IfNotSimd = 0> operator ComplexD () const { return(TensorRemove(_internal)); }; + // template = 0,IfNotSimd = 0> operator RealD () const { return(real(TensorRemove(_internal))); } + template = 0,IfNotSimd = 0> operator RealD () const { return TensorRemove(_internal); } + template = 0,IfNotSimd = 0> operator Integer () const { return Integer(TensorRemove(_internal)); } // convert from a something to a scalar via constructor of something arg template::value, T>::type* = nullptr > strong_inline iScalar operator = (T arg) diff --git a/lib/tensors/Tensor_peek.h b/lib/tensors/Tensor_peek.h index 5c8d4a93..9c037581 100644 --- a/lib/tensors/Tensor_peek.h +++ b/lib/tensors/Tensor_peek.h @@ -11,7 +11,7 @@ namespace Grid { //template inline ComplexD peekIndex(const ComplexD arg) { return arg;} //template inline RealF peekIndex(const RealF arg) { return arg;} //template inline RealD peekIndex(const RealD arg) { return arg;} - +#if 0 // Scalar peek, no indices template::TensorLevel == Level >::type * =nullptr> inline auto peekIndex(const iScalar &arg) -> iScalar @@ -88,6 +88,7 @@ template::T } return ret; } + // matrix template::TensorLevel != Level >::type * =nullptr> inline auto peekIndex(const iMatrix &arg) -> iMatrix(arg._internal[0][0])),N> @@ -119,6 +120,7 @@ template::T }} return ret; } +#endif } diff --git a/lib/tensors/Tensor_poke.h b/lib/tensors/Tensor_poke.h index a015478e..8a3fd69d 100644 --- a/lib/tensors/Tensor_poke.h +++ b/lib/tensors/Tensor_poke.h @@ -5,7 +5,7 @@ namespace Grid { ////////////////////////////////////////////////////////////////////////////// // Poke a specific index; ////////////////////////////////////////////////////////////////////////////// - +#if 0 // Scalar poke template::TensorLevel == Level >::type * =nullptr> inline void pokeIndex(iScalar &ret, const iScalar &arg) @@ -18,7 +18,7 @@ template::Te { ret._internal[i] = arg._internal; } -// Vector poke, two indices +//Matrix poke, two indices template::TensorLevel == Level >::type * =nullptr> inline void pokeIndex(iMatrix &ret, const iScalar &arg,int i,int j) { @@ -31,7 +31,6 @@ template::Te // scalar template::TensorLevel != Level >::type * =nullptr> inline void pokeIndex(iScalar &ret, const iScalar(ret._internal))> &arg) - { pokeIndex(ret._internal,arg._internal); } @@ -95,7 +94,7 @@ template::Te pokeIndex(ret._internal[ii][jj],arg._internal[ii][jj],i,j); }} } - +#endif } #endif diff --git a/lib/tensors/Tensor_trace.h b/lib/tensors/Tensor_trace.h index 749b292e..bcde670c 100644 --- a/lib/tensors/Tensor_trace.h +++ b/lib/tensors/Tensor_trace.h @@ -32,58 +32,327 @@ inline auto trace(const iScalar &arg) -> iScalar::TensorLevel != Level >::type * =nullptr> inline auto -traceIndex (const iScalar &arg) -> iScalar(arg._internal))> -{ - iScalar(arg._internal))> ret; - ret._internal=traceIndex(arg._internal); - return ret; -} -template::TensorLevel == Level >::type * =nullptr> inline auto -traceIndex (const iScalar &arg) -> iScalar -{ - return arg; -} -// If we hit the right index, return scalar and trace it with no further recursion -template::TensorLevel == Level >::type * =nullptr> inline -auto traceIndex(const iMatrix &arg) -> iScalar -{ - iScalar ret; - zeroit(ret._internal); - for(int i=0;i::TensorLevel != Level >::type * =nullptr> inline -auto traceIndex(const iMatrix &arg) -> iMatrix(arg._internal[0][0])),N> -{ - iMatrix(arg._internal[0][0])),N> ret; - for(int i=0;i(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::TensorLevel != Level >::type * =nullptr> inline -auto traceIndex(const iVector &arg) -> iVector(arg._internal[0])),N> -{ - iVector(arg._internal[0])),N> ret; - for(int i=0;i(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 +class TensorIndexRecursion { + + public: + //////////////////////////////////////////// + // Recursion for tracing a specific index + //////////////////////////////////////////// + template + static auto traceIndex(const iScalar arg) -> iScalar::traceIndex(arg._internal))> + { + iScalar::traceIndex(arg._internal))> ret; + ret._internal = TensorIndexRecursion::traceIndex(arg._internal); + return ret; } + template + static auto traceIndex(const iVector arg) -> iVector::traceIndex(arg._internal[0])),N> + { + iVector::traceIndex(arg._internal[0])),N> ret; + for(int i=0;i::traceIndex(arg._internal[i]); + } + return ret; + } + template + static auto traceIndex(const iMatrix arg) -> iMatrix::traceIndex(arg._internal[0][0])),N> + { + iMatrix::traceIndex(arg._internal[0][0])),N> ret; + for(int i=0;i::traceIndex(arg._internal[i][j]); + }} + return ret; + } + //////////////////////////////////////////// + // Recursion for peeking a specific index + //////////////////////////////////////////// + template + static auto peekIndex(const iScalar arg,int i) -> iScalar::peekIndex(arg._internal,0))> + { + iScalar::peekIndex(arg._internal,0))> ret; + ret._internal = TensorIndexRecursion::peekIndex(arg._internal,i); + return ret; + } + template + static auto peekIndex(const iScalar arg,int i,int j) -> iScalar::peekIndex(arg._internal,0,0))> + { + iScalar::peekIndex(arg._internal,0,0))> ret; + ret._internal = TensorIndexRecursion::peekIndex(arg._internal,i,j); + return ret; + } + + template + static auto peekIndex(const iVector arg,int ii) -> iVector::peekIndex(arg._internal[0],0)),N> + { + iVector::peekIndex(arg._internal[0],0)),N> ret; + for(int i=0;i::peekIndex(arg._internal[i],ii); + } + return ret; + } + template + static auto peekIndex(const iVector arg,int ii,int jj) -> iVector::peekIndex(arg._internal[0],0,0)),N> + { + iVector::peekIndex(arg._internal[0],0,0)),N> ret; + for(int i=0;i::peekIndex(arg._internal[i],ii,jj); + } + return ret; + } + + template + static auto peekIndex(const iMatrix arg,int ii) -> iMatrix::peekIndex(arg._internal[0][0],0)),N> + { + iMatrix::peekIndex(arg._internal[0][0],0)),N> ret; + for(int i=0;i::peekIndex(arg._internal[i][j],ii); + }} + return ret; + } + template + static auto peekIndex(const iMatrix arg,int ii,int jj) -> iMatrix::peekIndex(arg._internal[0][0],0,0)),N> + { + iMatrix::peekIndex(arg._internal[0][0],0,0)),N> ret; + for(int i=0;i::peekIndex(arg._internal[i][j],ii,jj); + }} + return ret; + } + //////////////////////////////////////////// + // Recursion for poking a specific index + //////////////////////////////////////////// + + template inline static + void pokeIndex(iScalar &ret, const iScalar::peekIndex(ret._internal,0))> &arg, int i) + { + TensorIndexRecursion::pokeIndex(ret._internal,arg._internal,i); + } + template inline static + void pokeIndex(iScalar &ret, const iScalar::peekIndex(ret._internal,0,0))> &arg, int i,int j) + { + TensorIndexRecursion::pokeIndex(ret._internal,arg._internal,i,j); + } + + template inline static + void pokeIndex(iVector &ret, const iVector::peekIndex(ret._internal,0)),N> &arg, int i) + { + for(int ii=0;ii::pokeIndex(ret._internal[ii],arg._internal[ii],i); + } + } + template inline static + void pokeIndex(iVector &ret, const iVector::peekIndex(ret._internal,0)),N> &arg, int i,int j) + { + for(int ii=0;ii::pokeIndex(ret._internal[ii],arg._internal[ii],i,j); + } + } + + template inline static + void pokeIndex(iMatrix &ret, const iMatrix::peekIndex(ret._internal,0)),N> &arg, int i) + { + for(int ii=0;ii::pokeIndex(ret._internal[ii][jj],arg._internal[ii][jj],i); + }} + } + template inline static + void pokeIndex(iMatrix &ret, const iMatrix::peekIndex(ret._internal,0)),N> &arg, int i,int j) + { + for(int ii=0;ii::pokeIndex(ret._internal[ii][jj],arg._internal[ii][jj],i,j); + }} + } + + //////////////////////////////////////////// + // Recursion for transposing a specific index + //////////////////////////////////////////// + template + static auto transposeIndex(const iScalar arg) -> iScalar + { + iScalar ret; + ret._internal = TensorIndexRecursion::transposeIndex(arg._internal); + return ret; + } + template + static auto transposeIndex(const iVector arg) -> iVector + { + iVector ret; + for(int i=0;i::transposeIndex(arg._internal[i]); + } + return ret; + } + template + static auto transposeIndex(const iMatrix arg) -> iMatrix + { + iMatrix ret; + for(int i=0;i::transposeIndex(arg._internal[i][j]); + }} + return ret; + } +}; + +//////////////////////////// +// strip const & ref quali's +//////////////////////////// +#define RemoveCRV(a) typename std::remove_const::type>::type +template<> +class TensorIndexRecursion<0> { + public: + + ///////////////////////////////////////// + // Ends recursion for trace (scalar/vector/matrix) + ///////////////////////////////////////// + template + static auto traceIndex(const iScalar arg) -> iScalar + { + iScalar ret; + ret._internal = arg._internal; + return ret; + } + template + static auto traceIndex(const iVector arg) -> iScalar + { + iScalar ret; + ret._internal=zero; + for(int i=0;i + static auto traceIndex(const iMatrix arg) -> iScalar + { + iScalar ret; + ret=zero; + for(int i=0;i + static auto transposeIndex(const iScalar arg) -> iScalar + { + iScalar ret; + ret._internal = arg._internal; + return ret; + } + template + static auto transposeIndex(const iVector arg) -> iVector + { + iScalar ret; + ret._internal=zero; + for(int i=0;i + static auto transposeIndex(const iMatrix arg) -> iMatrix + { + iScalar ret; + ret=zero; + for(int i=0;i + static auto peekIndex(const iVector arg,int ii) -> iScalar + { + iScalar ret; + ret._internal = arg._internal[ii]; + return ret; + } + template + static auto peekIndex(const iMatrix arg,int ii,int jj) -> iScalar + { + iScalar ret; + ret._internal = arg._internal[ii][jj]; + return ret; + } + // Vector poke, one index + template inline static + void pokeIndex(iVector &ret, const iScalar &arg,int i) + { + ret._internal[i] = arg._internal; + } + // Matrix poke two indices + template inline static + void pokeIndex(iMatrix &ret, const iScalar &arg,int i,int j) + { + ret._internal[i][j] = arg._internal; + } + +}; + +//////////////////////////////////////////////////////////////////////////////////////////////////////// +// External wrappers +//////////////////////////////////////////////////////////////////////////////////////////////////////// + +template inline auto traceIndex (const vtype &arg) -> RemoveCRV(TensorIndexRecursion::traceIndex(arg)) +{ + RemoveCRV(TensorIndexRecursion::traceIndex(arg)) ret; + ret=TensorIndexRecursion::traceIndex(arg); return ret; } +template inline auto transposeIndex (const vtype &arg) -> RemoveCRV(TensorIndexRecursion::transposeIndex(arg)) +{ + RemoveCRV(TensorIndexRecursion::transposeIndex(arg)) ret; + ret=TensorIndexRecursion::transposeIndex(arg); + return ret; +} +template inline auto peekIndex (const vtype &arg,int i) -> RemoveCRV(TensorIndexRecursion::peekIndex(arg,0)) +{ + RemoveCRV(TensorIndexRecursion::peekIndex(arg,0)) ret; + ret=TensorIndexRecursion::peekIndex(arg,i); + return ret; +} +template inline auto peekIndex (const vtype &arg,int i,int j) -> RemoveCRV(TensorIndexRecursion::peekIndex(arg,0,0)) +{ + RemoveCRV(TensorIndexRecursion::peekIndex(arg,0,0)) ret; + ret=TensorIndexRecursion::peekIndex(arg,i,j); + return ret; +} + +template inline +void pokeIndex (vtype &ret,const decltype(TensorIndexRecursion::peekIndex(ret,0)) &arg,int i) +{ + TensorIndexRecursion::pokeIndex(ret,arg,i); +} + +template inline +void pokeIndex (vtype &ret,const decltype(TensorIndexRecursion::peekIndex(ret,0,0)) &arg,int i,int j) +{ + TensorIndexRecursion::pokeIndex(ret,arg,i,j); +} + + +#undef RemoveCRV } #endif diff --git a/tests/Test_dwf_even_odd.cc b/tests/Test_dwf_even_odd.cc index 82b39f68..5ec73212 100644 --- a/tests/Test_dwf_even_odd.cc +++ b/tests/Test_dwf_even_odd.cc @@ -53,7 +53,7 @@ int main (int argc, char ** argv) random(RNG4,U[nn]); if ( nn>0 ) U[nn]=zero; - pokeIndex(Umu,U[nn],nn); + PokeIndex(Umu,U[nn],nn); } RealD mass=0.1; diff --git a/tests/Test_main.cc b/tests/Test_main.cc index 6bef8629..85450047 100644 --- a/tests/Test_main.cc +++ b/tests/Test_main.cc @@ -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(cVec,scalar,1); scalar=transpose(scalar); diff --git a/tests/Test_quenched_update.cc b/tests/Test_quenched_update.cc index 887a27a8..cb58144b 100644 --- a/tests/Test_quenched_update.cc +++ b/tests/Test_quenched_update.cc @@ -72,7 +72,7 @@ int main (int argc, char ** argv) } - pokeIndex(Umu,link,mu); + PokeIndex(Umu,link,mu); //reunitarise link; } diff --git a/tests/Test_wilson_even_odd.cc b/tests/Test_wilson_even_odd.cc index 16019c9c..91fbb90f 100644 --- a/tests/Test_wilson_even_odd.cc +++ b/tests/Test_wilson_even_odd.cc @@ -55,7 +55,7 @@ int main (int argc, char ** argv) Umu=zero; for(int nn=0;nn(Umu,U[nn],nn); + PokeIndex(Umu,U[nn],nn); } RealD mass=0.1;