portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-04 19:25:56 +01:00
Code Releases Activity

add benchmark to see whether matrix multiplication is slower than read from object

Browse Source
This commit is contained in:
david clarke 2023-06-27 21:28:26 -06:00
parent a7eabaad56
commit 9015c229dc
3 changed files with 188 additions and 164 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

149
Grid/qcd/smearing/HISQSmearing.h
View File

@ -2,9 +2,9 @@
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/smearing/StoutSmearing.h Source file: ./lib/qcd/smearing/HISQSmearing.h
Copyright (C) 2019 Copyright (C) 2023
Author: D. A. Clarke <clarke.davida@gmail.com> Author: D. A. Clarke <clarke.davida@gmail.com>
@ -56,7 +56,7 @@ template<class vobj> void gpermute(vobj & inout,int perm) {
/*! @brief signals that you want to go backwards in direction dir */ /*! @brief signals that you want to go backwards in direction dir */
inline int Back(const int dir) { inline int Back(const int dir) {
// generalShift will use BACKWARD_CONST to determine whether we step forward or // generalShift will use BACKWARD_CONST to determine whether we step forward or
// backward. Should work as long as BACKWARD_CONST > Nd. // backward. Should work as long as BACKWARD_CONST > Nd. Trick inspired by SIMULATeQCD.
return dir + BACKWARD_CONST; return dir + BACKWARD_CONST;
} }
@ -191,45 +191,33 @@ public:
for(int mu=0;mu<Nd;mu++) for(int mu=0;mu<Nd;mu++)
for(int nu=0;nu<Nd;nu++) { for(int nu=0;nu<Nd;nu++) {
if(mu==nu) continue; if(mu==nu) continue;
// x+mu auto SE0 = gStencil.GetEntry(s+0,ss); int x_p_mu = SE0->_offset;
// m+nu auto SE1 = gStencil.GetEntry(s+1,ss); int x_p_nu = SE1->_offset;
// x auto SE2 = gStencil.GetEntry(s+2,ss); int x = SE2->_offset;
// shift_munu; shift_munu[mu]= 1; shift_munu[nu]=-1; auto SE3 = gStencil.GetEntry(s+3,ss); int x_p_mu_m_nu = SE3->_offset;
// x-nu auto SE4 = gStencil.GetEntry(s+4,ss); int x_m_nu = SE4->_offset;
auto SE0 = gStencil.GetEntry(s+0,ss);
auto SE1 = gStencil.GetEntry(s+1,ss);
auto SE2 = gStencil.GetEntry(s+2,ss);
auto SE3 = gStencil.GetEntry(s+3,ss);
auto SE4 = gStencil.GetEntry(s+4,ss);
// Each offset corresponds to a site around the plaquette.
int o0 = SE0->_offset;
int o1 = SE1->_offset;
int o2 = SE2->_offset;
int o3 = SE3->_offset;
int o4 = SE4->_offset;
// When you're deciding whether to take an adjoint, the question is: how is the // When you're deciding whether to take an adjoint, the question is: how is the
// stored link oriented compared to the one you want? If I imagine myself travelling // stored link oriented compared to the one you want? If I imagine myself travelling
// with the to-be-updated link, I have two possible, alternative 3-link paths I can // with the to-be-updated link, I have two possible, alternative 3-link paths I can
// take, one starting by going to the left, the other starting by going to the right. // take, one starting by going to the left, the other starting by going to the right.
auto U0 = adj(U_v[o0](nu)); auto U0 = adj(U_v[x_p_mu](nu));
auto U1 = U_v[o1](mu); auto U1 = U_v[x_p_nu](mu) ;
auto U2 = U_v[o2](nu); auto U2 = U_v[x ](nu) ;
gpermute(U0,SE0->_permute); gpermute(U0,SE0->_permute);
gpermute(U1,SE1->_permute); gpermute(U1,SE1->_permute);
gpermute(U2,SE2->_permute); gpermute(U2,SE2->_permute);
auto U3 = U_v[o3](nu); auto U3 = U_v[x_p_mu_m_nu](nu) ;
auto U4 = U_v[o4](mu); auto U4 = U_v[x_m_nu ](mu) ;
auto U5 = adj(U_v[o4](nu)); auto U5 = adj(U_v[x_m_nu ](nu));
gpermute(U3,SE3->_permute); gpermute(U3,SE3->_permute);
gpermute(U4,SE4->_permute); gpermute(U4,SE4->_permute);
gpermute(U4,SE4->_permute);
// "left" "right" // "left" "right"
auto W = U2*U1*U0 + U5*U4*U3; auto W = U2*U1*U0 + U5*U4*U3;
@ -265,111 +253,8 @@ public:
~Smear_HISQ_Naik() {} ~Smear_HISQ_Naik() {}
void smear(LGF& u_smr, const LGF& U) const { // void smear(LGF& u_smr, const LGF& U) const {
// };
int depth = 1;
PaddedCell Ghost(depth,this->_grid);
LGF Ughost = Ghost.Exchange(u_smr);
GridBase *GhostGrid = Ughost.Grid();
LatticeComplex gplaq(GhostGrid);
LGF Ughost_naik(Ughost.Grid());
std::vector<Coordinate> shifts;
for(int mu=0;mu<Nd;mu++){
for(int nu=mu+1;nu<Nd;nu++){
// forward shifts
Coordinate x(Nd,0);
Coordinate shift_mu(Nd,0); shift_mu[mu]=1;
Coordinate shift_nu(Nd,0); shift_nu[nu]=1;
// push_back creates an element at the end of shifts and
// assigns the data in the argument to it.
shifts.push_back(shift_mu);
shifts.push_back(shift_nu);
shifts.push_back(x);
// reverse shifts
shift_nu[nu]=-1;
Coordinate shift_munu(Nd,0); shift_munu[mu]=1; shift_munu[nu]=-1;
shifts.push_back(shift_munu);
shifts.push_back(shift_nu); // in principle you don't need both of these grid points,
shifts.push_back(shift_nu); // but it helps the reader keep track of offsets
}
}
GeneralLocalStencil gStencil(GhostGrid,shifts);
Ughost_naik=Zero();
// Create the accessors, here U_v and U_fat_v
autoView(U_v , Ughost , CpuRead);
autoView(U_naik_v, Ughost_naik, CpuWrite);
// This is a loop over local sites.
for(int ss=0;ss<U_v.size();ss++){
// This is the stencil index. It increases as we make our way through the spacetime sites,
// plaquette orientations, and as we travel around a plaquette.
int s=0;
for(int mu=0;mu<Nd;mu++){
for(int nu=mu+1;nu<Nd;nu++){
// shift_mu; shift_mu[mu]=1
// shift_nu; shift_nu[nu]=1
// x
// shift_munu; shift_munu[mu]= 1; shift_munu[nu]=-1;
// shift_nu ; shift_nu[nu]=-1;
// shift_nu ; shift_nu[nu]=-1;
auto SE0 = gStencil.GetEntry(s+0,ss);
auto SE1 = gStencil.GetEntry(s+1,ss);
auto SE2 = gStencil.GetEntry(s+2,ss);
auto SE3 = gStencil.GetEntry(s+3,ss);
auto SE4 = gStencil.GetEntry(s+4,ss);
auto SE5 = gStencil.GetEntry(s+5,ss);
// Each offset corresponds to a site around the plaquette.
int o0 = SE0->_offset;
int o1 = SE1->_offset;
int o2 = SE2->_offset;
int o3 = SE3->_offset;
int o4 = SE4->_offset;
int o5 = SE5->_offset;
auto U0 = U_v[o0](nu);
auto U1 = adj(U_v[o1](mu));
auto U2 = adj(U_v[o2](nu));
gpermute(U0,SE0->_permute);
gpermute(U1,SE1->_permute);
gpermute(U2,SE2->_permute);
auto U3 = adj(U_v[o3](nu));
auto U4 = adj(U_v[o4](mu));
auto U5 = U_v[o5](nu);
gpermute(U3,SE3->_permute);
gpermute(U4,SE4->_permute);
gpermute(U5,SE5->_permute);
// Forward contribution from this orientation
auto W = U0*U1*U2;
U_naik_v[ss](mu) = U_naik_v[ss](mu) + W;
// Backward contribution from this orientation
W = U3*U4*U5;
U_naik_v[ss](mu) = U_naik_v[ss](mu) + W;
s=s+6;
}
}
}
// Here is my understanding of this part: The padded cell has its own periodic BCs, so
// if I take a step to the right at the right-most side of the cell, I end up on the
// left-most side. This means that the plaquettes in the padding are wrong. Luckily
// all we care about are the plaquettes in the cell, which we obtain from Extract.
u_smr = Ghost.Extract(Ughost_naik);
};
// void derivative(const GaugeField& Gauge) const { // void derivative(const GaugeField& Gauge) const {
// }; // };

125
benchmarks/Benchmark_su3mult_vs_lookup.cc Normal file
View File

@ -0,0 +1,125 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./benchmarks/Benchmark_su3mult_vs_lookup.cc
Copyright (C) 2023
Author: D. A. Clarke <clarke.davida@gmail.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/*
@file Benchmark_su3mult_vs_lookup.cc
@brief check to see whether su3 multiplication or lookup tables is faster
*/
#include <Grid/Grid.h>
using namespace Grid;
/*! @brief make the logger work like python print */
template<typename... Args>
inline std::string sjoin(Args&&... args) noexcept {
std::ostringstream msg;
(msg << ... << args);
return msg.str();
}
template <typename... Args>
inline void Grid_log(Args&&... args) {
std::string msg = sjoin(std::forward<Args>(args)...);
std::cout << GridLogMessage << msg << std::endl;
}
/*! @brief parameter file to easily adjust Nloop */
struct ConfParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(
ConfParameters,
int, Nloop);
template <class ReaderClass>
ConfParameters(Reader<ReaderClass>& Reader){
read(Reader, "parameters", *this);
}
};
int main (int argc, char** argv) {
// Params for the test.
int Ns = 8;
int Nt = 4;
int threads = GridThread::GetThreads();
std::string conf_in = "nersc.l8t4b3360";
Coordinate latt_size(Nd,0); latt_size[0]=Ns; latt_size[1]=Ns; latt_size[2]=Ns; latt_size[3]=Nt;
Grid_init(&argc,&argv);
Coordinate simd_layout = GridDefaultSimd(Nd,vComplexD::Nsimd());
Coordinate mpi_layout = GridDefaultMpi();
GridCartesian GRID(latt_size,simd_layout,mpi_layout);
Grid_log(" mpi = ",mpi_layout);
Grid_log(" simd = ",simd_layout);
Grid_log(" latt = ",latt_size);
Grid_log("threads = ",threads);
XmlReader Reader("mult_vs_lookup.xml",false, "grid");
ConfParameters param(Reader);
Grid_log(" Nloop = ",param.Nloop);
// Gauge field and accessor
LatticeGaugeField Umu(&GRID);
autoView(U_v, Umu, CpuRead);
// Read the configuration into Umu
FieldMetaData header;
NerscIO::readConfiguration(Umu, header, conf_in);
// Read in lattice sequentially, Nloop times
double lookupTime = 0.;
for(int i=0;i<param.Nloop;i++) {
double start = usecond();
for(int ss=0;ss<U_v.size();ss++)
for(int mu=0;mu<Nd;mu++) {
auto U1 = U_v[ss](mu);
}
double stop = usecond();
lookupTime += stop-start; // microseconds
}
Grid_log("Time to lookup: ",lookupTime,"[ms]");
// Raise a matrix to the power nmat, for each link.
auto U1 = U_v[0](0);
for(int nmat=1;nmat<8;nmat++) {
double multTime = 0.;
for(int i=0;i<param.Nloop;i++) {
double start=usecond();
for(int ss=0;ss<U_v.size();ss++)
for(int mu=0;mu<Nd;mu++) {
auto U2 = U1;
for(int j=1;j<nmat;j++) {
U2 *= U1;
}
}
double stop=usecond();
multTime += stop-start;
}
Grid_log("Time to multiply ",nmat," matrices: ",lookupTime," [ms]");
}
Grid_finalize();
}

78
tests/smearing/Test_fatLinks.cc
View File

@ -1,11 +1,35 @@
/* /*************************************************************************************
* Test_fatLinks.cc
* Grid physics library, www.github.com/paboyle/Grid
* D. Clarke
* Source file: ./tests/smearing/Test_fatLinks.cc
* Test the various constructs used to make fat links.
* Copyright (C) 2023
*/
Author: D. A. Clarke <clarke.davida@gmail.com>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/*
@file Test_fatLinks.cc
@brief test of the HISQ smearing
*/
#include <Grid/Grid.h> #include <Grid/Grid.h>
#include <Grid/lattice/PaddedCell.h> #include <Grid/lattice/PaddedCell.h>
@ -14,7 +38,7 @@
using namespace Grid; using namespace Grid;
// Make the logger work like Python print() /*! @brief make the logger work like python print */
template<typename... Args> template<typename... Args>
inline std::string sjoin(Args&&... args) noexcept { inline std::string sjoin(Args&&... args) noexcept {
std::ostringstream msg; std::ostringstream msg;
@ -27,32 +51,26 @@ inline void Grid_log(Args&&... args) {
std::cout << GridLogMessage << msg << std::endl; std::cout << GridLogMessage << msg << std::endl;
} }
struct fatParams: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(
fatParams,
std::string, conf_in,
std::string, conf_out);
template <class ReaderClass>
fatParams(Reader<ReaderClass>& Reader){
read(Reader, "parameters", *this);
}
};
// //
// one method: input --> fat // one method: input --> fat
// another : input --> long (naik) // another : input --> long (naik)
// another : input --> unitarize // another : input --> unitarize
// //
int main (int argc, char **argv) int main (int argc, char** argv) {
{
// Params for the test.
int Ns = 8;
int Nt = 4;
Coordinate latt_size(Nd,0); latt_size[0]=Ns; latt_size[1]=Ns; latt_size[2]=Ns; latt_size[3]=Nt;
std::string conf_in = "nersc.l8t4b3360";
std::string conf_out = "nersc.l8t4b3360.3link";
// Initialize the Grid // Initialize the Grid
Grid_init(&argc,&argv); Grid_init(&argc,&argv);
Coordinate latt_size = GridDefaultLatt();
Coordinate simd_layout = GridDefaultSimd(Nd,vComplexD::Nsimd()); Coordinate simd_layout = GridDefaultSimd(Nd,vComplexD::Nsimd());
Coordinate mpi_layout = GridDefaultMpi(); Coordinate mpi_layout = GridDefaultMpi();
Grid_log("mpi = ",mpi_layout); Grid_log(" mpi = ",mpi_layout);
Grid_log("simd = ",simd_layout); Grid_log("simd = ",simd_layout);
Grid_log("latt = ",latt_size); Grid_log("latt = ",latt_size);
GridCartesian GRID(latt_size,simd_layout,mpi_layout); GridCartesian GRID(latt_size,simd_layout,mpi_layout);
@ -61,19 +79,15 @@ int main (int argc, char **argv)
LatticeGaugeField Umu(&GRID); LatticeGaugeField Umu(&GRID);
LatticeGaugeField U_smr(&GRID); LatticeGaugeField U_smr(&GRID);
// Read in the parameter file
XmlReader Reader("fatParams.xml",false, "grid");
fatParams param(Reader);
FieldMetaData header;
// Read the configuration into Umu // Read the configuration into Umu
NerscIO::readConfiguration(Umu, header, param.conf_in); FieldMetaData header;
NerscIO::readConfiguration(Umu, header, conf_in);
// Smear Umu and store result in U_smr // Smear Umu and store result in U_smr
Smear_HISQ_fat<LatticeGaugeField> hisq_fat(&GRID,1/8.,0.,1/16.,1/64.,1/384.,0.); Smear_HISQ_fat<LatticeGaugeField> hisq_fat(&GRID,1/8.,0.,1/16.,1/64.,1/384.,0.);
hisq_fat.smear(U_smr,Umu); hisq_fat.smear(U_smr,Umu);
NerscIO::writeConfiguration(U_smr,param.conf_out,"HISQ"); NerscIO::writeConfiguration(U_smr,conf_out,"HISQ");
// Test a C-style instantiation // Test a C-style instantiation
double path_coeff[6] = {1, 2, 3, 4, 5, 6}; double path_coeff[6] = {1, 2, 3, 4, 5, 6};