diff --git a/tests/Test_stencil.cc b/tests/Test_stencil.cc index cc49ab22..5949fa63 100644 --- a/tests/Test_stencil.cc +++ b/tests/Test_stencil.cc @@ -8,6 +8,10 @@ int main (int argc, char ** argv) { Grid_init(&argc,&argv); + // typedef LatticeColourMatrix Field; + typedef LatticeComplex Field; + typedef typename Field::vector_object vobj; + typedef typename vobj::scalar_object sobj; std::vector latt_size = GridDefaultLatt(); std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); @@ -18,23 +22,40 @@ int main (int argc, char ** argv) GridCartesian Fine(latt_size,simd_layout,mpi_layout); GridRedBlackCartesian rbFine(latt_size,simd_layout,mpi_layout); GridParallelRNG fRNG(&Fine); + // fRNG.SeedRandomDevice(); std::vector seeds({1,2,3,4}); fRNG.SeedFixedIntegers(seeds); - LatticeColourMatrix Foo(&Fine); - LatticeColourMatrix Bar(&Fine); - LatticeColourMatrix Check(&Fine); - LatticeColourMatrix Diff(&Fine); - + Field Foo(&Fine); + Field Bar(&Fine); + Field Check(&Fine); + Field Diff(&Fine); + LatticeComplex lex(&Fine); + + lex = zero; random(fRNG,Foo); gaussian(fRNG,Bar); + /* + Integer stride =1000; + { + double nrm; + LatticeComplex coor(&Fine); + + for(int d=0;d directions(npoint,dir); @@ -48,8 +69,8 @@ int main (int argc, char ** argv) ocoor[dir]=(ocoor[dir]+disp)%Fine._rdimensions[dir]; } - std::vector > comm_buf(myStencil._unified_buffer_size); - SimpleCompressor compress; + std::vector > comm_buf(myStencil._unified_buffer_size); + SimpleCompressor compress; myStencil.HaloExchange(Foo,comm_buf,compress); Bar = Cshift(Foo,dir,disp); @@ -75,9 +96,114 @@ int main (int argc, char ** argv) Real nrm = norm2(Diff); std::cout< coor(4); + for(coor[3]=0;coor[3] 0){ + std::cout <<"Coor (" << coor[0]<<","< directions(npoint,dir); + std::vector displacements(npoint,disp); + + CartesianStencil EStencil(&rbFine,npoint,Even,directions,displacements); + CartesianStencil OStencil(&rbFine,npoint,Odd,directions,displacements); + + std::vector ocoor(4); + for(int o=0;o > Ecomm_buf(EStencil._unified_buffer_size); + std::vector > Ocomm_buf(OStencil._unified_buffer_size); + + SimpleCompressor compress; + + EStencil.HaloExchange(EFoo,Ecomm_buf,compress); + OStencil.HaloExchange(OFoo,Ocomm_buf,compress); + + Bar = Cshift(Foo,dir,disp); + + if ( disp & 0x1 ) { + ECheck.checkerboard = Even; + OCheck.checkerboard = Odd; + } else { + ECheck.checkerboard = Odd; + OCheck.checkerboard = Even; + } + // Implement a stencil code that should agree with that darn cshift! + for(int i=0;ioSites();i++){ + int permute_type; + StencilEntry *SE; + SE = EStencil.GetEntry(permute_type,0,i); + std::cout << "Even source "<< i<<" -> " <_offset << " "<< SE->_is_local<_is_local && SE->_permute ) + permute(OCheck._odata[i],EFoo._odata[SE->_offset],permute_type); + else if (SE->_is_local) + OCheck._odata[i] = EFoo._odata[SE->_offset]; + else + OCheck._odata[i] = Ecomm_buf[SE->_offset]; + } + for(int i=0;ioSites();i++){ + int permute_type; + StencilEntry *SE; + SE = OStencil.GetEntry(permute_type,0,i); + std::cout << "ODD source "<< i<<" -> " <_offset << " "<< SE->_is_local<_is_local && SE->_permute ) + permute(ECheck._odata[i],OFoo._odata[SE->_offset],permute_type); + else if (SE->_is_local) + ECheck._odata[i] = OFoo._odata[SE->_offset]; + else + ECheck._odata[i] = Ocomm_buf[SE->_offset]; + } + + setCheckerboard(Check,ECheck); + setCheckerboard(Check,OCheck); + + Real nrmC = norm2(Check); + Real nrmB = norm2(Bar); + Diff = Check-Bar; + Real nrm = norm2(Diff); + std::cout< coor(4); for(coor[3]=0;coor[3] 0){ - printf("Coor (%d %d %d %d) \t rc %d%d \t %le (%le,%le) %le\n", - coor[0],coor[1],coor[2],coor[3],r,c, - nn, - real(check()()(r,c)), - imag(check()()(r,c)), - real(bar()()(r,c)) - ); - } - snrmC=snrmC+real(conjugate(check()()(r,c))*check()()(r,c)); - snrmB=snrmB+real(conjugate(bar()()(r,c))*bar()()(r,c)); - snrm=snrm+nn; - }} + sobj ddiff; + ddiff = check -bar; + diff =norm2(ddiff); + if ( diff > 0){ + std::cout <<"Coor (" << coor[0]<<","<