diff --git a/benchmarks/Benchmark_dwf.cc b/benchmarks/Benchmark_dwf.cc index 10e4521b..76625973 100644 --- a/benchmarks/Benchmark_dwf.cc +++ b/benchmarks/Benchmark_dwf.cc @@ -163,6 +163,24 @@ int main (int argc, char ** argv) Dw.Report(); } + if (1) { // Naive wilson dag implementation + ref = zero; + for (int mu = 0; mu < Nd; mu++) { + // ref = src - Gamma(Gamma::GammaX)* src ; // 1+gamma_x + tmp = U[mu] * Cshift(src, mu + 1, 1); + for (int i = 0; i < ref._odata.size(); i++) { + ref._odata[i] += tmp._odata[i] + Gamma(Gmu[mu]) * tmp._odata[i]; + } + + tmp = adj(U[mu]) * src; + tmp = Cshift(tmp, mu + 1, -1); + for (int i = 0; i < ref._odata.size(); i++) { + ref._odata[i] += tmp._odata[i] - Gamma(Gmu[mu]) * tmp._odata[i]; + } + } + ref = -0.5 * ref; + } + if (1) { @@ -245,6 +263,33 @@ int main (int argc, char ** argv) std::cout<::Dhop to naive wilson implementation Dag to verify correctness" << std::endl; + sDw.Dhop(ssrc,sresult,1); + sum=0; + for(int x=0;x site({s,x,y,z,t}); + SpinColourVector normal, simd; + peekSite(normal,ref,site); + peekSite(simd,sresult,site); + sum=sum+norm2(normal-simd); + if (norm2(normal-simd) > 1.0e-6 ) { + std::cout << "site "<1.0e-4) { setCheckerboard(ssrc,ssrc_o); setCheckerboard(ssrc,ssrc_e); std::cout<< ssrc << std::endl; } + + // Check the dag + std::cout << GridLogMessage << "Compare WilsonFermion5D::DhopEO to Dhop to verify correctness" << std::endl; + pickCheckerboard(Even,ssrc_e,ssrc); + pickCheckerboard(Odd,ssrc_o,ssrc); + sDw.DhopEO(ssrc_o,sr_e,DaggerYes); + sDw.DhopOE(ssrc_e,sr_o,DaggerYes); + sDw.Dhop (ssrc ,sresult,DaggerYes); + + pickCheckerboard(Even,ssrc_e,sresult); + pickCheckerboard(Odd ,ssrc_o,sresult); + ssrc_e = ssrc_e - sr_e; + error = norm2(ssrc_e); + + std::cout<1.0e-4) { + setCheckerboard(ssrc,ssrc_o); + setCheckerboard(ssrc,ssrc_e); + std::cout<< ssrc << std::endl; + } + + } } - if (1) - { // Naive wilson dag implementation - ref = zero; - for(int mu=0;mu { std::cout << GridLogDebug << "a = "<< a << " b_pred = "<< b_pred << " b = "<< b << std::endl; std::cout << GridLogDebug << "qq = "<< qq << " d = "<< d << " c = "<< c << std::endl; - // Stopping condition if (cp <= rsq) { SolverTimer.Stop(); diff --git a/lib/simd/Grid_avx.h b/lib/simd/Grid_avx.h index 36360102..77927b10 100644 --- a/lib/simd/Grid_avx.h +++ b/lib/simd/Grid_avx.h @@ -514,7 +514,7 @@ namespace Optimization { template static inline __m256 tRotate(__m256 in){ __m256 tmp = Permute::Permute0(in); - __m256 ret; + __m256 ret = in; if ( n > 3 ) { _mm256_alignr_epi32_grid(ret,in,tmp,n); } else { @@ -526,7 +526,7 @@ namespace Optimization { template static inline __m256d tRotate(__m256d in){ __m256d tmp = Permute::Permute0(in); - __m256d ret; + __m256d ret = in; if ( n > 1 ) { _mm256_alignr_epi64_grid(ret,in,tmp,n); } else {