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

Merge remote-tracking branch 'upstream/master'

Browse Source
This commit is contained in:
neo 2015-05-29 11:41:39 +09:00
commit 575e6001f3
11 changed files with 335 additions and 363 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
INSTALL
View File

@ -1 +1 @@
/usr/share/automake-1.14/INSTALL
/opt/local/share/automake-1.15/INSTALL

31
benchmarks/Grid_memory_bandwidth.cc
View File

@ -27,8 +27,8 @@ int main (int argc, char ** argv)
for(int lat=4;lat<=32;lat+=4){
std::vector<int> latt_size ({lat,lat,lat,lat});
std::vector<int> latt_size ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
int vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
//GridParallelRNG pRNG(&Grid); pRNG.SeedRandomDevice();
@ -47,8 +47,8 @@ int main (int argc, char ** argv)
double stop=usecond();
double time = (stop-start)/Nloop*1000;
double flops=lat*lat*lat*lat*Nvec*2;// mul,add
double bytes=3*lat*lat*lat*lat*Nvec*sizeof(Real);
double flops=vol*Nvec*2;// mul,add
double bytes=3*vol*Nvec*sizeof(Real);
std::cout<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t"<<flops/time<<std::endl;
}
@ -61,8 +61,8 @@ int main (int argc, char ** argv)
for(int lat=4;lat<=32;lat+=4){
std::vector<int> latt_size ({lat,lat,lat,lat});
std::vector<int> latt_size ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
int vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
//GridParallelRNG pRNG(&Grid); pRNG.SeedRandomDevice();
@ -79,8 +79,8 @@ int main (int argc, char ** argv)
double stop=usecond();
double time = (stop-start)/Nloop*1000;
double flops=lat*lat*lat*lat*Nvec*2;// mul,add
double bytes=3*lat*lat*lat*lat*Nvec*sizeof(Real);
double flops=vol*Nvec*2;// mul,add
double bytes=3*vol*Nvec*sizeof(Real);
std::cout<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t"<<flops/time<<std::endl;
}
@ -92,7 +92,8 @@ int main (int argc, char ** argv)
for(int lat=4;lat<=32;lat+=4){
std::vector<int> latt_size ({lat,lat,lat,lat});
std::vector<int> latt_size ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
int vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
@ -111,8 +112,8 @@ int main (int argc, char ** argv)
double stop=usecond();
double time = (stop-start)/Nloop*1000;
double bytes=2*lat*lat*lat*lat*Nvec*sizeof(Real);
double flops=lat*lat*lat*lat*Nvec*1;// mul
double bytes=2*vol*Nvec*sizeof(Real);
double flops=vol*Nvec*1;// mul
std::cout <<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t"<<flops/time<<std::endl;
}
@ -125,8 +126,8 @@ int main (int argc, char ** argv)
for(int lat=4;lat<=32;lat+=4){
std::vector<int> latt_size ({lat,lat,lat,lat});
std::vector<int> latt_size ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
int vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
//GridParallelRNG pRNG(&Grid); pRNG.SeedRandomDevice();
@ -144,8 +145,8 @@ int main (int argc, char ** argv)
double stop=usecond();
double time = (stop-start)/Nloop*1000;
double bytes=lat*lat*lat*lat*Nvec*sizeof(Real);
double flops=lat*lat*lat*lat*Nvec*2;// mul,add
double bytes=vol*Nvec*sizeof(Real);
double flops=vol*Nvec*2;// mul,add
std::cout<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t"<<flops/time<<std::endl;
}

31
benchmarks/Grid_su3.cc
View File

@ -24,8 +24,8 @@ int main (int argc, char ** argv)
for(int lat=2;lat<=24;lat+=2){
std::vector<int> latt_size ({lat,lat,lat,lat});
std::vector<int> latt_size ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
int vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
// GridParallelRNG pRNG(&Grid); pRNG.SeedRandomDevice();
@ -40,9 +40,9 @@ int main (int argc, char ** argv)
double stop=usecond();
double time = (stop-start)/Nloop*1000.0;
double bytes=3.0*lat*lat*lat*lat*Nc*Nc*sizeof(Complex);
double footprint=2.0*lat*lat*lat*lat*Nc*Nc*sizeof(Complex);
double flops=Nc*Nc*(6.0+8.0+8.0)*lat*lat*lat*lat;
double bytes=3.0*vol*Nc*Nc*sizeof(Complex);
double footprint=2.0*vol*Nc*Nc*sizeof(Complex);
double flops=Nc*Nc*(6.0+8.0+8.0)*vol;
std::cout<<std::setprecision(3) << lat<<"\t\t"<<footprint<<" \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
}
@ -56,7 +56,8 @@ int main (int argc, char ** argv)
for(int lat=2;lat<=24;lat+=2){
std::vector<int> latt_size ({lat,lat,lat,lat});
std::vector<int> latt_size ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
int vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
// GridParallelRNG pRNG(&Grid); pRNG.SeedRandomDevice();
@ -72,8 +73,8 @@ int main (int argc, char ** argv)
double stop=usecond();
double time = (stop-start)/Nloop*1000.0;
double bytes=3*lat*lat*lat*lat*Nc*Nc*sizeof(Complex);
double flops=Nc*Nc*(6+8+8)*lat*lat*lat*lat;
double bytes=3*vol*Nc*Nc*sizeof(Complex);
double flops=Nc*Nc*(6+8+8)*vol;
std::cout<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
}
@ -86,7 +87,8 @@ int main (int argc, char ** argv)
for(int lat=2;lat<=24;lat+=2){
std::vector<int> latt_size ({lat,lat,lat,lat});
std::vector<int> latt_size ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
int vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
// GridParallelRNG pRNG(&Grid); pRNG.SeedRandomDevice();
@ -102,8 +104,8 @@ int main (int argc, char ** argv)
double stop=usecond();
double time = (stop-start)/Nloop*1000.0;
double bytes=3*lat*lat*lat*lat*Nc*Nc*sizeof(Complex);
double flops=Nc*Nc*(6+8+8)*lat*lat*lat*lat;
double bytes=3*vol*Nc*Nc*sizeof(Complex);
double flops=Nc*Nc*(6+8+8)*vol;
std::cout<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
}
@ -116,7 +118,8 @@ int main (int argc, char ** argv)
for(int lat=2;lat<=24;lat+=2){
std::vector<int> latt_size ({lat,lat,lat,lat});
std::vector<int> latt_size ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
int vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
// GridParallelRNG pRNG(&Grid); pRNG.SeedRandomDevice();
@ -132,8 +135,8 @@ int main (int argc, char ** argv)
double stop=usecond();
double time = (stop-start)/Nloop*1000.0;
double bytes=3*lat*lat*lat*lat*Nc*Nc*sizeof(Complex);
double flops=Nc*Nc*(8+8+8)*lat*lat*lat*lat;
double bytes=3*vol*Nc*Nc*sizeof(Complex);
double flops=Nc*Nc*(8+8+8)*vol;
std::cout<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
}

13
configure vendored
View File

@ -2574,7 +2574,7 @@ test -n "$target_alias" &&
NONENONEs,x,x, &&
program_prefix=${target_alias}-
am__api_version='1.14'
am__api_version='1.15'
# Find a good install program. We prefer a C program (faster),
# so one script is as good as another. But avoid the broken or
@ -2746,8 +2746,8 @@ test "$program_suffix" != NONE &&
ac_script='s/[\\$]/&&/g;s/;s,x,x,$//'
program_transform_name=`$as_echo "$program_transform_name" | sed "$ac_script"`
# expand $ac_aux_dir to an absolute path
am_aux_dir=`cd $ac_aux_dir && pwd`
# Expand $ac_aux_dir to an absolute path.
am_aux_dir=`cd "$ac_aux_dir" && pwd`
if test x"${MISSING+set}" != xset; then
case $am_aux_dir in
@ -2766,7 +2766,7 @@ else
$as_echo "$as_me: WARNING: 'missing' script is too old or missing" >&2;}
fi
if test x"${install_sh}" != xset; then
if test x"${install_sh+set}" != xset; then
case $am_aux_dir in
*\ * | *\ *)
install_sh="\${SHELL} '$am_aux_dir/install-sh'" ;;
@ -3094,8 +3094,8 @@ MAKEINFO=${MAKEINFO-"${am_missing_run}makeinfo"}
# <http://lists.gnu.org/archive/html/automake/2012-07/msg00014.html>
mkdir_p='$(MKDIR_P)'
# We need awk for the "check" target. The system "awk" is bad on
# some platforms.
# We need awk for the "check" target (and possibly the TAP driver). The
# system "awk" is bad on some platforms.
# Always define AMTAR for backward compatibility. Yes, it's still used
# in the wild :-( We should find a proper way to deprecate it ...
AMTAR='$${TAR-tar}'
@ -3154,6 +3154,7 @@ END
fi
ac_config_headers="$ac_config_headers lib/Grid_config.h"
# Check whether --enable-silent-rules was given.

1
lib/Grid_init.cc
View File

@ -54,7 +54,6 @@ namespace Grid {
const std::vector<int> &GridDefaultLatt(void) {return Grid_default_latt;};
const std::vector<int> &GridDefaultMpi(void) {return Grid_default_mpi;};
////////////////////////////////////////////////////////////
// Command line parsing assist for stock controls
////////////////////////////////////////////////////////////

2
lib/Makefile.am
View File

@ -96,7 +96,7 @@ nobase_include_HEADERS = algorithms/approx/bigfloat.h \
simd/Grid_vector_types.h \
simd/Grid_sse4.h \
simd/Grid_avx.h \
simd/Grid_knc.h
simd/Grid_avx512.h

43
lib/cshift/Grid_cshift_common.h
View File

@ -27,9 +27,11 @@ Gather_plane_simple (const Lattice<vobj> &rhs,std::vector<cobj,alignedAllocator<
int so = plane*rhs._grid->_ostride[dimension]; // base offset for start of plane
int e1=rhs._grid->_slice_nblock[dimension];
int e2=rhs._grid->_slice_block[dimension];
PARALLEL_NESTED_LOOP2
for(int n=0;n<rhs._grid->_slice_nblock[dimension];n++){
for(int b=0;b<rhs._grid->_slice_block[dimension];b++){
for(int n=0;n<e1;n++){
for(int b=0;b<e2;b++){
int o = n*rhs._grid->_slice_stride[dimension];
int bo = n*rhs._grid->_slice_block[dimension];
int ocb=1<<rhs._grid->CheckerBoardFromOindex(o+b);// Could easily be a table lookup
@ -54,10 +56,12 @@ Gather_plane_extract(const Lattice<vobj> &rhs,std::vector<typename cobj::scalar_
}
int so = plane*rhs._grid->_ostride[dimension]; // base offset for start of plane
int e1=rhs._grid->_slice_nblock[dimension];
int e2=rhs._grid->_slice_block[dimension];
PARALLEL_NESTED_LOOP2
for(int n=0;n<rhs._grid->_slice_nblock[dimension];n++){
for(int b=0;b<rhs._grid->_slice_block[dimension];b++){
for(int n=0;n<e1;n++){
for(int b=0;b<e2;b++){
int o=n*rhs._grid->_slice_stride[dimension];
int offset = b+n*rhs._grid->_slice_block[dimension];
@ -103,9 +107,11 @@ template<class vobj> void Scatter_plane_simple (Lattice<vobj> &rhs,std::vector<v
int so = plane*rhs._grid->_ostride[dimension]; // base offset for start of plane
int e1=rhs._grid->_slice_nblock[dimension];
int e2=rhs._grid->_slice_block[dimension];
PARALLEL_NESTED_LOOP2
for(int n=0;n<rhs._grid->_slice_nblock[dimension];n++){
for(int b=0;b<rhs._grid->_slice_block[dimension];b++){
for(int n=0;n<e1;n++){
for(int b=0;b<e2;b++){
int o =n*rhs._grid->_slice_stride[dimension];
int bo =n*rhs._grid->_slice_block[dimension];
int ocb=1<<rhs._grid->CheckerBoardFromOindex(o+b);// Could easily be a table lookup
@ -129,10 +135,11 @@ PARALLEL_NESTED_LOOP2
int so = plane*rhs._grid->_ostride[dimension]; // base offset for start of plane
int e1=rhs._grid->_slice_nblock[dimension];
int e2=rhs._grid->_slice_block[dimension];
PARALLEL_NESTED_LOOP2
for(int n=0;n<rhs._grid->_slice_nblock[dimension];n++){
for(int b=0;b<rhs._grid->_slice_block[dimension];b++){
for(int n=0;n<e1;n++){
for(int b=0;b<e2;b++){
int o = n*rhs._grid->_slice_stride[dimension];
int offset = b+n*rhs._grid->_slice_block[dimension];
int ocb=1<<rhs._grid->CheckerBoardFromOindex(o+b);
@ -156,10 +163,12 @@ template<class vobj> void Copy_plane(Lattice<vobj>& lhs,Lattice<vobj> &rhs, int
int ro = rplane*rhs._grid->_ostride[dimension]; // base offset for start of plane
int lo = lplane*lhs._grid->_ostride[dimension]; // base offset for start of plane
int e1=rhs._grid->_slice_nblock[dimension]; // clearly loop invariant for icpc
int e2=rhs._grid->_slice_block[dimension];
PARALLEL_NESTED_LOOP2
for(int n=0;n<rhs._grid->_slice_nblock[dimension];n++){
for(int b=0;b<rhs._grid->_slice_block[dimension];b++){
for(int n=0;n<e1;n++){
for(int b=0;b<e2;b++){
int o =n*rhs._grid->_slice_stride[dimension]+b;
int ocb=1<<lhs._grid->CheckerBoardFromOindex(o);
@ -185,10 +194,12 @@ template<class vobj> void Copy_plane_permute(Lattice<vobj>& lhs,Lattice<vobj> &r
int ro = rplane*rhs._grid->_ostride[dimension]; // base offset for start of plane
int lo = lplane*lhs._grid->_ostride[dimension]; // base offset for start of plane
int e1=rhs._grid->_slice_nblock[dimension];
int e2=rhs._grid->_slice_block [dimension];
PARALLEL_NESTED_LOOP2
for(int n=0;n<rhs._grid->_slice_nblock[dimension];n++){
for(int b=0;b<rhs._grid->_slice_block [dimension];b++){
for(int n=0;n<e1;n++){
for(int b=0;b<e2;b++){
int o =n*rhs._grid->_slice_stride[dimension];
int ocb=1<<lhs._grid->CheckerBoardFromOindex(o+b);

75
lib/simd/Grid_avx.h
View File

@ -14,7 +14,23 @@
#endif
namespace Optimization {
template<class vtype>
union uconv {
__m256 f;
vtype v;
};
union u256f {
__m256 v;
float f[8];
};
union u256d {
__m256d v;
double f[4];
};
struct Vsplat{
//Complex float
inline __m256 operator()(float a, float b){
@ -54,7 +70,6 @@ namespace Optimization {
};
struct Vstream{
//Float
inline void operator()(float * a, __m256 b){
@ -68,8 +83,6 @@ namespace Optimization {
};
struct Vset{
// Complex float
inline __m256 operator()(Grid::ComplexF *a){
@ -92,7 +105,6 @@ namespace Optimization {
return _mm256_set_epi32(a[7],a[6],a[5],a[4],a[3],a[2],a[1],a[0]);
}
};
template <typename Out_type, typename In_type>
@ -106,9 +118,6 @@ namespace Optimization {
}
};
/////////////////////////////////////////////////////
// Arithmetic operations
/////////////////////////////////////////////////////
@ -170,7 +179,6 @@ namespace Optimization {
}
};
struct MultComplex{
// Complex float
inline __m256 operator()(__m256 a, __m256 b){
@ -207,7 +215,6 @@ namespace Optimization {
IF IMM0[3] = 0
THEN DEST[255:192]=SRC2[191:128] ELSE DEST[255:192]=SRC2[255:192] FI; // Ox5 r<->i ; 0xC unchanged
*/
__m256d ymm0,ymm1,ymm2;
ymm0 = _mm256_shuffle_pd(a,a,0x0); // ymm0 <- ar ar, ar,ar b'00,00
ymm0 = _mm256_mul_pd(ymm0,b); // ymm0 <- ar bi, ar br
@ -247,7 +254,6 @@ namespace Optimization {
}
};
struct Conj{
// Complex single
inline __m256 operator()(__m256 in){
@ -292,18 +298,13 @@ namespace Optimization {
}
};
//////////////////////////////////////////////
// Some Template specialization
//////////////////////////////////////////////
template < typename vtype >
void permute(vtype &a, vtype &b, int perm) {
union {
__m256 f;
vtype v;
} conv;
void permute(vtype &a,vtype b, int perm) {
uconv<vtype> conv;
conv.v = b;
switch (perm){
// 8x32 bits=>3 permutes
@ -313,24 +314,20 @@ namespace Optimization {
default: assert(0); break;
}
a = conv.v;
}
//Complex float Reduce
template<>
inline Grid::ComplexF Reduce<Grid::ComplexF, __m256>::operator()(__m256 in){
__m256 v1,v2;
union {
__m256 v;
float f[8];
} conv;
Optimization::permute(v1,in,0); // sse 128; paired complex single
v1 = _mm256_add_ps(v1,in);
Optimization::permute(v2,v1,1); // avx 256; quad complex single
v1 = _mm256_add_ps(v1,v2);
conv.v = v1;
u256f conv; conv.v = v1;
return Grid::ComplexF(conv.f[0],conv.f[1]);
}
//Real float Reduce
template<>
inline Grid::RealF Reduce<Grid::RealF, __m256>::operator()(__m256 in){
@ -341,7 +338,8 @@ namespace Optimization {
v1 = _mm256_add_ps(v1,v2);
Optimization::permute(v2,v1,2);
v1 = _mm256_add_ps(v1,v2);
return v1[0];
u256f conv; conv.v=v1;
return conv.f[0];
}
@ -351,7 +349,8 @@ namespace Optimization {
__m256d v1;
Optimization::permute(v1,in,0); // sse 128; paired complex single
v1 = _mm256_add_pd(v1,in);
return Grid::ComplexD(v1[0],v1[1]);
u256d conv; conv.v = v1;
return Grid::ComplexD(conv.f[0],conv.f[1]);
}
//Real double Reduce
@ -362,7 +361,8 @@ namespace Optimization {
v1 = _mm256_add_pd(v1,in);
Optimization::permute(v2,v1,1);
v1 = _mm256_add_pd(v1,v2);
return v1[0];
u256d conv; conv.v = v1;
return conv.f[0];
}
//Integer Reduce
@ -390,22 +390,9 @@ namespace Grid {
_mm_prefetch(ptr+i+512,_MM_HINT_T0);
}
}
template < typename VectorSIMD >
inline void Gpermute(VectorSIMD &y,const VectorSIMD &b, int perm ) {
union {
__m256 f;
decltype(VectorSIMD::v) v;
} conv;
conv.v = b.v;
switch(perm){
case 3: break; //empty for AVX1/2
case 2: conv.f = _mm256_shuffle_ps(conv.f,conv.f,_MM_SHUFFLE(2,3,0,1)); break;
case 1: conv.f = _mm256_shuffle_ps(conv.f,conv.f,_MM_SHUFFLE(1,0,3,2)); break;
case 0: conv.f = _mm256_permute2f128_ps(conv.f,conv.f,0x01); break;
default: assert(0); break;
}
y.v=conv.v;
Optimization::permute(y.v,b.v,perm);
};
// Function name aliases

0
lib/simd/Grid_knc.h → lib/simd/Grid_avx512.h
View File

109
lib/simd/Grid_sse4.h
View File

@ -10,6 +10,21 @@
#include <pmmintrin.h>
namespace Optimization {
template<class vtype>
union uconv {
__m128 f;
vtype v;
};
union u128f {
__m128 v;
float f[4];
};
union u128d {
__m128d v;
double f[2];
};
struct Vsplat{
//Complex float
@ -50,7 +65,6 @@ namespace Optimization {
};
struct Vstream{
//Float
inline void operator()(float * a, __m128 b){
@ -64,8 +78,6 @@ namespace Optimization {
};
struct Vset{
// Complex float
inline __m128 operator()(Grid::ComplexF *a){
@ -102,9 +114,6 @@ namespace Optimization {
}
};
/////////////////////////////////////////////////////
// Arithmetic operations
/////////////////////////////////////////////////////
@ -138,7 +147,6 @@ namespace Optimization {
}
};
struct MultComplex{
// Complex float
inline __m128 operator()(__m128 a, __m128 b){
@ -177,7 +185,6 @@ namespace Optimization {
}
};
struct Conj{
// Complex single
inline __m128 operator()(__m128 in){
@ -216,57 +223,61 @@ namespace Optimization {
__m128d tmp = _mm_shuffle_pd(in,in,0x1);
return _mm_addsub_pd(_mm_setzero_pd(),tmp); // r,-i
}
};
//////////////////////////////////////////////
// Some Template specialization
template < typename vtype >
void permute(vtype &a, vtype b, int perm) {
uconv<vtype> conv;
conv.v = b;
switch(perm){
case 3: break; //empty for SSE4
case 2: break; //empty for SSE4
case 1: conv.f = _mm_shuffle_ps(conv.f,conv.f,_MM_SHUFFLE(2,3,0,1)); break;
case 0: conv.f = _mm_shuffle_ps(conv.f,conv.f,_MM_SHUFFLE(1,0,3,2)); break;
default: assert(0); break;
}
a=conv.v;
};
//Complex float Reduce
template<>
inline Grid::ComplexF Reduce<Grid::ComplexF, __m128>::operator()(__m128 in){
union {
__m128 v1;
float f[4];
} u128;
u128.v1 = _mm_add_ps(in, _mm_shuffle_ps(in,in, 0b01001110)); // FIXME Prefer to use _MM_SHUFFLE macros
return Grid::ComplexF(u128.f[0], u128.f[1]);
__m128 v1; // two complex
Optimization::permute(v1,in,0);
v1= _mm_add_ps(v1,in);
u128f conv; conv.v=v1;
return Grid::ComplexF(conv.f[0],conv.f[1]);
}
//Real float Reduce
template<>
inline Grid::RealF Reduce<Grid::RealF, __m128>::operator()(__m128 in){
// FIXME Hack
const Grid::RealF * ptr = (const Grid::RealF *) &in;
Grid::RealF ret = 0;
for(int i=0;i< 4 ;i++){ // 4 number of simd lanes for float
ret = ret+ptr[i];
}
return ret;
__m128 v1,v2; // quad single
Optimization::permute(v1,in,0);
v1= _mm_add_ps(v1,in);
Optimization::permute(v2,v1,1);
v1 = _mm_add_ps(v1,v2);
u128f conv; conv.v=v1;
return conv.f[0];
}
//Complex double Reduce
template<>
inline Grid::ComplexD Reduce<Grid::ComplexD, __m128d>::operator()(__m128d in){
printf("Reduce : Missing good complex double implementation -> FIX\n");
return Grid::ComplexD(in[0], in[1]); // inefficient
u128d conv; conv.v = in;
return Grid::ComplexD(conv.f[0],conv.f[1]);
}
//Real double Reduce
template<>
inline Grid::RealD Reduce<Grid::RealD, __m128d>::operator()(__m128d in){
// FIXME Hack
const Grid::RealD * ptr =(const Grid::RealD *) &in;
Grid::RealD ret = 0;
for(int i=0;i< 2 ;i++){// 2 number of simd lanes for float
ret = ret+ptr[i];
}
return ret;
__m128d v1;
Optimization::permute(v1,in,0); // avx 256; quad double
v1 = _mm_add_pd(v1,in);
u128d conv; conv.v = v1;
return conv.f[0];
}
//Integer Reduce
@ -276,12 +287,6 @@ namespace Optimization {
printf("Reduce : Missing integer implementation -> FIX\n");
assert(0);
}
}
//////////////////////////////////////////////////////////////////////////////////////
@ -292,27 +297,13 @@ namespace Grid {
typedef __m128d SIMD_Dtype; // Double precision type
typedef __m128i SIMD_Itype; // Integer type
inline void v_prefetch0(int size, const char *ptr){}; // prefetch utilities
// Gpermute function
template < typename VectorSIMD >
inline void Gpermute(VectorSIMD &y,const VectorSIMD &b, int perm ) {
union {
__m128 f;
decltype(VectorSIMD::v) v;
} conv;
conv.v = b.v;
switch(perm){
case 3: break; //empty for SSE4
case 2: break; //empty for SSE4
case 1: conv.f = _mm_shuffle_ps(conv.f,conv.f,_MM_SHUFFLE(2,3,0,1)); break;
case 0: conv.f = _mm_shuffle_ps(conv.f,conv.f,_MM_SHUFFLE(1,0,3,2)); break;
default: assert(0); break;
}
y.v=conv.v;
};
Optimization::permute(y.v,b.v,perm);
}
// Function name aliases

391
lib/simd/Grid_vector_types.h
View File

@ -14,7 +14,7 @@
#include "Grid_avx.h"
#endif
#if defined AVX512
#include "Grid_knc.h"
#include "Grid_avx512.h"
#endif
#if defined QPX
#include "Grid_qpx.h"
@ -32,13 +32,24 @@ namespace Grid {
// type alias used to simplify the syntax of std::enable_if
template <typename T> using Invoke = typename T::type;
template <typename Condition, typename ReturnType> using EnableIf = Invoke<std::enable_if<Condition::value, ReturnType>>;
template <typename Condition, typename ReturnType> using NotEnableIf= Invoke<std::enable_if<!Condition::value, ReturnType>>;
template <typename Condition, typename ReturnType> using EnableIf = Invoke<std::enable_if<Condition::value, ReturnType> >;
template <typename Condition, typename ReturnType> using NotEnableIf= Invoke<std::enable_if<!Condition::value, ReturnType> >;
////////////////////////////////////////////////////////
// Check for complexity with type traits
template <typename T> struct is_complex : std::false_type {};
template < typename T > struct is_complex< std::complex<T> >: std::true_type {};
template <typename T> struct is_complex : public std::false_type {};
template <> struct is_complex<std::complex<double> >: public std::true_type {};
template <> struct is_complex<std::complex<float> > : public std::true_type {};
template <typename T> using IfReal = Invoke<std::enable_if<std::is_floating_point<T>::value,int> > ;
template <typename T> using IfComplex = Invoke<std::enable_if<is_complex<T>::value,int> > ;
template <typename T> using IfInteger = Invoke<std::enable_if<std::is_integral<T>::value,int> > ;
template <typename T> using IfNotReal = Invoke<std::enable_if<!std::is_floating_point<T>::value,int> > ;
template <typename T> using IfNotComplex = Invoke<std::enable_if<!is_complex<T>::value,int> > ;
template <typename T> using IfNotInteger = Invoke<std::enable_if<!std::is_integral<T>::value,int> > ;
////////////////////////////////////////////////////////
// Define the operation templates functors
// general forms to allow for vsplat syntax
@ -54,11 +65,9 @@ namespace Grid {
Out unary(Input src, Operation op){
return op(src);
}
///////////////////////////////////////////////
/*
@brief Grid_simd class for the SIMD vector type operations
@ -73,27 +82,28 @@ namespace Grid {
Vector_type v;
static inline int Nsimd(void) { return sizeof(Vector_type)/sizeof(Scalar_type);}
// Constructors
Grid_simd & operator = ( Zero & z){
vzero(*this);
return (*this);
}
Grid_simd& operator=(const Grid_simd&& rhs){v=rhs.v;return *this;};
Grid_simd& operator=(const Grid_simd& rhs){v=rhs.v;return *this;}; //faster than not declaring it and leaving to the compiler
Grid_simd()=default;
Grid_simd(const Grid_simd& rhs):v(rhs.v){}; //compiles in movaps
Grid_simd(const Grid_simd& rhs) :v(rhs.v){}; //compiles in movaps
Grid_simd(const Grid_simd&& rhs):v(rhs.v){};
/////////////////////////////
// Constructors
/////////////////////////////
Grid_simd & operator = ( Zero & z){
vzero(*this);
return (*this);
}
//Enable if complex type
template < class S = Scalar_type >
template < typename S = Scalar_type >
Grid_simd(const typename std::enable_if< is_complex < S >::value, S>::type a){
vsplat(*this,a);
};
Grid_simd(const Real a){
vsplat(*this,Scalar_type(a));
@ -107,86 +117,16 @@ namespace Grid {
friend inline void sub (Grid_simd * __restrict__ y,const Grid_simd * __restrict__ l,const Grid_simd *__restrict__ r){ *y = (*l) - (*r); }
friend inline void add (Grid_simd * __restrict__ y,const Grid_simd * __restrict__ l,const Grid_simd *__restrict__ r){ *y = (*l) + (*r); }
friend inline void mac (Grid_simd *__restrict__ y,const Scalar_type *__restrict__ a,const Grid_simd *__restrict__ x){ *y = (*a)*(*x)+(*y); };
friend inline void mult(Grid_simd *__restrict__ y,const Scalar_type *__restrict__ l,const Grid_simd *__restrict__ r){ *y = (*l) * (*r); }
friend inline void sub (Grid_simd *__restrict__ y,const Scalar_type *__restrict__ l,const Grid_simd *__restrict__ r){ *y = (*l) - (*r); }
friend inline void add (Grid_simd *__restrict__ y,const Scalar_type *__restrict__ l,const Grid_simd *__restrict__ r){ *y = (*l) + (*r); }
friend inline void mac (Grid_simd *__restrict__ y,const Grid_simd *__restrict__ a,const Scalar_type *__restrict__ x){ *y = (*a)*(*x)+(*y); };
friend inline void mult(Grid_simd *__restrict__ y,const Grid_simd *__restrict__ l,const Scalar_type *__restrict__ r){ *y = (*l) * (*r); }
friend inline void sub (Grid_simd *__restrict__ y,const Grid_simd *__restrict__ l,const Scalar_type *__restrict__ r){ *y = (*l) - (*r); }
friend inline void add (Grid_simd *__restrict__ y,const Grid_simd *__restrict__ l,const Scalar_type *__restrict__ r){ *y = (*l) + (*r); }
//not for integer types...
template < class S = Scalar_type, NotEnableIf<std::is_integral < S >, int> = 0 >
friend inline Grid_simd adj(const Grid_simd &in){ return conjugate(in); }
///////////////////////////////////////////////
// Initialise to 1,0,i for the correct types
///////////////////////////////////////////////
// For complex types
template < class S = Scalar_type, EnableIf<is_complex < S >, int> = 0 >
friend inline void vone(Grid_simd &ret) { vsplat(ret,1.0,0.0); }
template < class S = Scalar_type, EnableIf<is_complex < S >, int> = 0 >
friend inline void vzero(Grid_simd &ret) { vsplat(ret,0.0,0.0); }// use xor?
template < class S = Scalar_type, EnableIf<is_complex < S >, int> = 0 >
friend inline void vcomplex_i(Grid_simd &ret){ vsplat(ret,0.0,1.0);}
// if not complex overload here
template < class S = Scalar_type, EnableIf<std::is_floating_point < S >,int> = 0 >
friend inline void vone(Grid_simd &ret) { vsplat(ret,1.0); }
template < class S = Scalar_type, EnableIf<std::is_floating_point < S >,int> = 0 >
friend inline void vzero(Grid_simd &ret) { vsplat(ret,0.0); }
// For integral types
template < class S = Scalar_type, EnableIf<std::is_integral < S >, int> = 0 >
friend inline void vone(Grid_simd &ret) { vsplat(ret,1); }
template < class S = Scalar_type, EnableIf<std::is_integral < S >, int> = 0 >
friend inline void vzero(Grid_simd &ret) { vsplat(ret,0); }
template < class S = Scalar_type, EnableIf<std::is_integral < S >, int> = 0 >
friend inline void vtrue (Grid_simd &ret){vsplat(ret,0xFFFFFFFF);}
template < class S = Scalar_type, EnableIf<std::is_integral < S >, int> = 0 >
friend inline void vfalse(Grid_simd &ret){vsplat(ret,0);}
////////////////////////////////////
// Arithmetic operator overloads +,-,*
////////////////////////////////////
friend inline Grid_simd operator + (Grid_simd a, Grid_simd b)
{
Grid_simd ret;
ret.v = binary<Vector_type>(a.v, b.v, SumSIMD());
return ret;
};
friend inline Grid_simd operator - (Grid_simd a, Grid_simd b)
{
Grid_simd ret;
ret.v = binary<Vector_type>(a.v, b.v, SubSIMD());
return ret;
};
// Distinguish between complex types and others
template < class S = Scalar_type, EnableIf<is_complex < S >, int> = 0 >
friend inline Grid_simd operator * (Grid_simd a, Grid_simd b)
{
Grid_simd ret;
ret.v = binary<Vector_type>(a.v,b.v, MultComplexSIMD());
return ret;
};
// Real/Integer types
template < class S = Scalar_type, NotEnableIf<is_complex < S >, int> = 0 >
friend inline Grid_simd operator * (Grid_simd a, Grid_simd b)
{
Grid_simd ret;
ret.v = binary<Vector_type>(a.v,b.v, MultSIMD());
return ret;
};
////////////////////////////////////////////////////////////////////////
// FIXME: gonna remove these load/store, get, set, prefetch
////////////////////////////////////////////////////////////////////////
@ -194,28 +134,6 @@ namespace Grid {
ret.v = unary<Vector_type>(a, VsetSIMD());
}
///////////////////////
// Splat
///////////////////////
// overload if complex
template < class S = Scalar_type >
friend inline void vsplat(Grid_simd &ret, EnableIf<is_complex < S >, S> c){
Real a = real(c);
Real b = imag(c);
vsplat(ret,a,b);
}
// this is only for the complex version
template < class S = Scalar_type, EnableIf<is_complex < S >, int> = 0 >
friend inline void vsplat(Grid_simd &ret,Real a, Real b){
ret.v = binary<Vector_type>(a, b, VsplatSIMD());
}
//if real fill with a, if complex fill with a in the real part (first function above)
friend inline void vsplat(Grid_simd &ret,Real a){
ret.v = unary<Vector_type>(a, VsplatSIMD());
}
///////////////////////
// Vstore
///////////////////////
@ -223,19 +141,6 @@ namespace Grid {
binary<void>(ret.v, (Real*)a, VstoreSIMD());
}
///////////////////////
// Vstream
///////////////////////
template < class S = Scalar_type, NotEnableIf<std::is_integral < S >, int> = 0 >
friend inline void vstream(Grid_simd &out,const Grid_simd &in){
binary<void>((Real*)&out.v, in.v, VstreamSIMD());
}
template < class S = Scalar_type, EnableIf<std::is_integral < S >, int> = 0 >
friend inline void vstream(Grid_simd &out,const Grid_simd &in){
out=in;
}
///////////////////////
// Vprefetch
///////////////////////
@ -244,7 +149,6 @@ namespace Grid {
_mm_prefetch((const char*)&v.v,_MM_HINT_T0);
}
///////////////////////
// Reduce
///////////////////////
@ -265,63 +169,6 @@ namespace Grid {
return a*b;
}
///////////////////////
// Conjugate
///////////////////////
template < class S = Scalar_type, EnableIf<is_complex < S >, int> = 0 >
friend inline Grid_simd conjugate(const Grid_simd &in){
Grid_simd ret ;
ret.v = unary<Vector_type>(in.v, ConjSIMD());
return ret;
}
template < class S = Scalar_type, NotEnableIf<is_complex < S >, int> = 0 >
friend inline Grid_simd conjugate(const Grid_simd &in){
return in; // for real objects
}
///////////////////////
// timesMinusI
///////////////////////
template < class S = Scalar_type, EnableIf<is_complex < S >, int> = 0 >
friend inline void timesMinusI( Grid_simd &ret,const Grid_simd &in){
ret.v = binary<Vector_type>(in.v, ret.v, TimesMinusISIMD());
}
template < class S = Scalar_type, EnableIf<is_complex < S >, int> = 0 >
friend inline Grid_simd timesMinusI(const Grid_simd &in){
Grid_simd ret;
timesMinusI(ret,in);
return ret;
}
template < class S = Scalar_type, NotEnableIf<is_complex < S >, int> = 0 >
friend inline Grid_simd timesMinusI(const Grid_simd &in){
return in;
}
///////////////////////
// timesI
///////////////////////
template < class S = Scalar_type, EnableIf<is_complex < S >, int> = 0 >
friend inline void timesI(Grid_simd &ret,const Grid_simd &in){
ret.v = binary<Vector_type>(in.v, ret.v, TimesISIMD());
}
template < class S = Scalar_type, EnableIf<is_complex < S >, int> = 0 >
friend inline Grid_simd timesI(const Grid_simd &in){
Grid_simd ret;
timesI(ret,in);
return ret;
}
template < class S = Scalar_type, NotEnableIf<is_complex < S >, int> = 0 >
friend inline Grid_simd timesI(const Grid_simd &in){
return in;
}
///////////////////////
// Unary negation
///////////////////////
@ -346,9 +193,6 @@ namespace Grid {
return *this;
}
////////////////////////////////////////////////////////////////////
// General permute; assumes vector length is same across
// all subtypes; may not be a good assumption, but could
@ -359,48 +203,183 @@ namespace Grid {
Gpermute<Grid_simd>(y,b,perm);
}
};// end of Grid_simd class definition
///////////////////////
// Splat
///////////////////////
// this is only for the complex version
template <class S, class V, IfComplex<S> =0, class ABtype>
inline void vsplat(Grid_simd<S,V> &ret,ABtype a, ABtype b){
ret.v = binary<V>(a, b, VsplatSIMD());
}
template<class scalar_type, class vector_type >
inline Grid_simd< scalar_type, vector_type> innerProduct(const Grid_simd< scalar_type, vector_type> & l, const Grid_simd< scalar_type, vector_type> & r)
// overload if complex
template <class S,class V> inline void vsplat(Grid_simd<S,V> &ret, EnableIf<is_complex < S >, S> c) {
Real a = real(c);
Real b = imag(c);
vsplat(ret,a,b);
}
//if real fill with a, if complex fill with a in the real part (first function above)
template <class S,class V>
inline void vsplat(Grid_simd<S,V> &ret,NotEnableIf<is_complex< S>,S> a){
ret.v = unary<V>(a, VsplatSIMD());
}
//////////////////////////
///////////////////////////////////////////////
// Initialise to 1,0,i for the correct types
///////////////////////////////////////////////
// For complex types
template <class S,class V, IfComplex<S> = 0 > inline void vone(Grid_simd<S,V> &ret) { vsplat(ret,S(1.0,0.0)); }
template <class S,class V, IfComplex<S> = 0 > inline void vzero(Grid_simd<S,V> &ret) { vsplat(ret,S(0.0,0.0)); }// use xor?
template <class S,class V, IfComplex<S> = 0 > inline void vcomplex_i(Grid_simd<S,V> &ret){ vsplat(ret,S(0.0,1.0));}
// if not complex overload here
template <class S,class V, IfReal<S> = 0 > inline void vone (Grid_simd<S,V> &ret){ vsplat(ret,1.0); }
template <class S,class V, IfReal<S> = 0 > inline void vzero(Grid_simd<S,V> &ret) { vsplat(ret,0.0); }
// For integral types
template <class S,class V,IfInteger<S> = 0 > inline void vone(Grid_simd<S,V> &ret) {vsplat(ret,1); }
template <class S,class V,IfInteger<S> = 0 > inline void vzero(Grid_simd<S,V> &ret) {vsplat(ret,0); }
template <class S,class V,IfInteger<S> = 0 > inline void vtrue (Grid_simd<S,V> &ret){vsplat(ret,0xFFFFFFFF);}
template <class S,class V,IfInteger<S> = 0 > inline void vfalse(Grid_simd<S,V> &ret){vsplat(ret,0);}
template<class S,class V> inline void zeroit(Grid_simd<S,V> &z){ vzero(z);}
///////////////////////
// Vstream
///////////////////////
template <class S,class V, IfNotInteger<S> = 0 >
inline void vstream(Grid_simd<S,V> &out,const Grid_simd<S,V> &in){
binary<void>((Real*)&out.v, in.v, VstreamSIMD());
}
template <class S,class V, IfInteger<S> = 0 >
inline void vstream(Grid_simd<S,V> &out,const Grid_simd<S,V> &in){
out=in;
}
////////////////////////////////////
// Arithmetic operator overloads +,-,*
////////////////////////////////////
template<class S,class V> inline Grid_simd<S,V> operator + (Grid_simd<S,V> a, Grid_simd<S,V> b) {
Grid_simd<S,V> ret;
ret.v = binary<V>(a.v, b.v, SumSIMD());
return ret;
};
template<class S,class V> inline Grid_simd<S,V> operator - (Grid_simd<S,V> a, Grid_simd<S,V> b) {
Grid_simd<S,V> ret;
ret.v = binary<V>(a.v, b.v, SubSIMD());
return ret;
};
// Distinguish between complex types and others
template<class S,class V, IfComplex<S> = 0 > inline Grid_simd<S,V> operator * (Grid_simd<S,V> a, Grid_simd<S,V> b) {
Grid_simd<S,V> ret;
ret.v = binary<V>(a.v,b.v, MultComplexSIMD());
return ret;
};
// Real/Integer types
template<class S,class V, IfNotComplex<S> = 0 > inline Grid_simd<S,V> operator * (Grid_simd<S,V> a, Grid_simd<S,V> b) {
Grid_simd<S,V> ret;
ret.v = binary<V>(a.v,b.v, MultSIMD());
return ret;
};
///////////////////////
// Conjugate
///////////////////////
template <class S,class V, IfComplex<S> = 0 >
inline Grid_simd<S,V> conjugate(const Grid_simd<S,V> &in){
Grid_simd<S,V> ret ;
ret.v = unary<V>(in.v, ConjSIMD());
return ret;
}
template <class S,class V, IfNotComplex<S> = 0 > inline Grid_simd<S,V> conjugate(const Grid_simd<S,V> &in){
return in; // for real objects
}
//Suppress adj for integer types... // odd; why conjugate above but not adj??
template < class S, class V, IfNotInteger<S> = 0 >
inline Grid_simd<S,V> adj(const Grid_simd<S,V> &in){ return conjugate(in); }
///////////////////////
// timesMinusI
///////////////////////
template<class S,class V,IfComplex<S> = 0 >
inline void timesMinusI( Grid_simd<S,V> &ret,const Grid_simd<S,V> &in){
ret.v = binary<V>(in.v, ret.v, TimesMinusISIMD());
}
template<class S,class V,IfComplex<S> = 0 >
inline Grid_simd<S,V> timesMinusI(const Grid_simd<S,V> &in){
Grid_simd<S,V> ret;
timesMinusI(ret,in);
return ret;
}
template<class S,class V,IfNotComplex<S> = 0 >
inline Grid_simd<S,V> timesMinusI(const Grid_simd<S,V> &in){
return in;
}
///////////////////////
// timesI
///////////////////////
template<class S,class V,IfComplex<S> = 0 >
inline void timesI(Grid_simd<S,V> &ret,const Grid_simd<S,V> &in){
ret.v = binary<V>(in.v, ret.v, TimesISIMD());
}
template<class S,class V,IfComplex<S> = 0 >
inline Grid_simd<S,V> timesI(const Grid_simd<S,V> &in){
Grid_simd<S,V> ret;
timesI(ret,in);
return ret;
}
template<class S,class V,IfNotComplex<S> = 0 >
inline Grid_simd<S,V> timesI(const Grid_simd<S,V> &in){
return in;
}
/////////////////////
// Inner, outer
/////////////////////
template<class S, class V >
inline Grid_simd< S, V> innerProduct(const Grid_simd< S, V> & l, const Grid_simd< S, V> & r)
{
return conjugate(l)*r;
}
template<class scalar_type, class vector_type >
inline void zeroit(Grid_simd< scalar_type, vector_type> &z){ vzero(z);}
template<class scalar_type, class vector_type >
inline Grid_simd< scalar_type, vector_type> outerProduct(const Grid_simd< scalar_type, vector_type> &l, const Grid_simd< scalar_type, vector_type>& r)
template<class S, class V >
inline Grid_simd< S, V> outerProduct(const Grid_simd< S, V> &l, const Grid_simd< S, V> & r)
{
return l*r;
}
template<class scalar_type, class vector_type >
inline Grid_simd< scalar_type, vector_type> trace(const Grid_simd< scalar_type, vector_type> &arg){
template<class S, class V >
inline Grid_simd< S, V> trace(const Grid_simd< S, V> &arg){
return arg;
}
///////////////////////////////
// Define available types
///////////////////////////////
typedef Grid_simd< float , SIMD_Ftype > vRealF;
typedef Grid_simd< double , SIMD_Dtype > vRealD;
typedef Grid_simd< std::complex< float > , SIMD_Ftype > vComplexF;
typedef Grid_simd< std::complex< double >, SIMD_Dtype > vComplexD;
typedef Grid_simd< Integer , SIMD_Itype > vInteger;
}
#endif