diff --git a/README b/README
index 41b66b6b..17e92fa0 100644
--- a/README
+++ b/README
@@ -33,6 +33,8 @@ MPI parallelism is UNIMPLEMENTED and for now only OpenMP and SIMD parallelism is
by setting variables in the command line or in the environment. Here
is are examples:
+ ./configure CXX=clang++ CXXFLAGS="-std=c++11 -O3 -msse4" --enable-simd=SSE4
+
./configure CXX=clang++ CXXFLAGS="-std=c++11 -O3 -mavx" --enable-simd=AVX1
./configure CXX=clang++ CXXFLAGS="-std=c++11 -O3 -mavx2" --enable-simd=AVX2
diff --git a/README.md b/README.md
index 80714d76..e18ca474 100644
--- a/README.md
+++ b/README.md
@@ -37,6 +37,8 @@ MPI, OpenMP, and SIMD parallelism are present in the library.
by setting variables in the command line or in the environment. Here
are examples:
+ ./configure CXX=clang++ CXXFLAGS="-std=c++11 -O3 -msse4" --enable-simd=SSE4
+
./configure CXX=clang++ CXXFLAGS="-std=c++11 -O3 -mavx" --enable-simd=AVX1
./configure CXX=clang++ CXXFLAGS="-std=c++11 -O3 -mavx2" --enable-simd=AVX2
diff --git a/configure b/configure
index b42143e2..6e27ab11 100755
--- a/configure
+++ b/configure
@@ -4503,6 +4503,11 @@ _ACEOF
# Checks for library functions.
+echo
+echo Checking libraries
+echo :::::::::::::::::::::::::::::::::::::::::::
+
+
for ac_func in gettimeofday
do :
ac_fn_cxx_check_func "$LINENO" "gettimeofday" "ac_cv_func_gettimeofday"
@@ -4515,6 +4520,105 @@ fi
done
+{ $as_echo "$as_me:${as_lineno-$LINENO}: checking for __gmpf_init in -lgmp" >&5
+$as_echo_n "checking for __gmpf_init in -lgmp... " >&6; }
+if ${ac_cv_lib_gmp___gmpf_init+:} false; then :
+ $as_echo_n "(cached) " >&6
+else
+ ac_check_lib_save_LIBS=$LIBS
+LIBS="-lgmp $LIBS"
+cat confdefs.h - <<_ACEOF >conftest.$ac_ext
+/* end confdefs.h. */
+
+/* Override any GCC internal prototype to avoid an error.
+ Use char because int might match the return type of a GCC
+ builtin and then its argument prototype would still apply. */
+#ifdef __cplusplus
+extern "C"
+#endif
+char __gmpf_init ();
+int
+main ()
+{
+return __gmpf_init ();
+ ;
+ return 0;
+}
+_ACEOF
+if ac_fn_cxx_try_link "$LINENO"; then :
+ ac_cv_lib_gmp___gmpf_init=yes
+else
+ ac_cv_lib_gmp___gmpf_init=no
+fi
+rm -f core conftest.err conftest.$ac_objext \
+ conftest$ac_exeext conftest.$ac_ext
+LIBS=$ac_check_lib_save_LIBS
+fi
+{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $ac_cv_lib_gmp___gmpf_init" >&5
+$as_echo "$ac_cv_lib_gmp___gmpf_init" >&6; }
+if test "x$ac_cv_lib_gmp___gmpf_init" = xyes; then :
+ cat >>confdefs.h <<_ACEOF
+#define HAVE_LIBGMP 1
+_ACEOF
+
+ LIBS="-lgmp $LIBS"
+
+else
+ as_fn_error $? "GNU Multiple Precision GMP library was not found in your system.
+Please install or provide the correct path to your installation
+Info at: http://www.gmplib.org" "$LINENO" 5
+fi
+
+
+{ $as_echo "$as_me:${as_lineno-$LINENO}: checking for mpfr_init in -lmpfr" >&5
+$as_echo_n "checking for mpfr_init in -lmpfr... " >&6; }
+if ${ac_cv_lib_mpfr_mpfr_init+:} false; then :
+ $as_echo_n "(cached) " >&6
+else
+ ac_check_lib_save_LIBS=$LIBS
+LIBS="-lmpfr $LIBS"
+cat confdefs.h - <<_ACEOF >conftest.$ac_ext
+/* end confdefs.h. */
+
+/* Override any GCC internal prototype to avoid an error.
+ Use char because int might match the return type of a GCC
+ builtin and then its argument prototype would still apply. */
+#ifdef __cplusplus
+extern "C"
+#endif
+char mpfr_init ();
+int
+main ()
+{
+return mpfr_init ();
+ ;
+ return 0;
+}
+_ACEOF
+if ac_fn_cxx_try_link "$LINENO"; then :
+ ac_cv_lib_mpfr_mpfr_init=yes
+else
+ ac_cv_lib_mpfr_mpfr_init=no
+fi
+rm -f core conftest.err conftest.$ac_objext \
+ conftest$ac_exeext conftest.$ac_ext
+LIBS=$ac_check_lib_save_LIBS
+fi
+{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $ac_cv_lib_mpfr_mpfr_init" >&5
+$as_echo "$ac_cv_lib_mpfr_mpfr_init" >&6; }
+if test "x$ac_cv_lib_mpfr_mpfr_init" = xyes; then :
+ cat >>confdefs.h <<_ACEOF
+#define HAVE_LIBMPFR 1
+_ACEOF
+
+ LIBS="-lmpfr $LIBS"
+
+else
+ as_fn_error $? "GNU Multiple Precision MPFR library was not found in your system.
+Please install or provide the correct path to your installation
+Info at: http://www.mpfr.org/" "$LINENO" 5
+fi
+
diff --git a/configure.ac b/configure.ac
index 14170f4e..93e2b574 100644
--- a/configure.ac
+++ b/configure.ac
@@ -3,7 +3,7 @@
#
# Project Grid package
#
-# Time-stamp: <2015-05-18 17:14:20 neo>
+# Time-stamp: <2015-05-19 13:51:08 neo>
AC_PREREQ([2.69])
AC_INIT([Grid], [1.0], [paboyle@ph.ed.ac.uk])
@@ -46,8 +46,22 @@ AC_TYPE_UINT32_T
AC_TYPE_UINT64_T
# Checks for library functions.
+echo
+echo Checking libraries
+echo :::::::::::::::::::::::::::::::::::::::::::
+
+
AC_CHECK_FUNCS([gettimeofday])
+AC_CHECK_LIB([gmp],[__gmpf_init],,
+ [AC_MSG_ERROR(GNU Multiple Precision GMP library was not found in your system.
+Please install or provide the correct path to your installation
+Info at: http://www.gmplib.org)])
+
+AC_CHECK_LIB([mpfr],[mpfr_init],,
+ [AC_MSG_ERROR(GNU Multiple Precision MPFR library was not found in your system.
+Please install or provide the correct path to your installation
+Info at: http://www.mpfr.org/)])
diff --git a/lib/Grid_config.h b/lib/Grid_config.h
index e1674850..78582f3e 100644
--- a/lib/Grid_config.h
+++ b/lib/Grid_config.h
@@ -36,6 +36,12 @@
/* Define to 1 if you have the <inttypes.h> header file. */
#define HAVE_INTTYPES_H 1
+/* Define to 1 if you have the `gmp' library (-lgmp). */
+#define HAVE_LIBGMP 1
+
+/* Define to 1 if you have the `mpfr' library (-lmpfr). */
+#define HAVE_LIBMPFR 1
+
/* Define to 1 if you have the <malloc.h> header file. */
#define HAVE_MALLOC_H 1
diff --git a/lib/Grid_config.h.in b/lib/Grid_config.h.in
index 0ce09cf5..b7f56d5b 100644
--- a/lib/Grid_config.h.in
+++ b/lib/Grid_config.h.in
@@ -35,6 +35,12 @@
/* Define to 1 if you have the <inttypes.h> header file. */
#undef HAVE_INTTYPES_H
+/* Define to 1 if you have the `gmp' library (-lgmp). */
+#undef HAVE_LIBGMP
+
+/* Define to 1 if you have the `mpfr' library (-lmpfr). */
+#undef HAVE_LIBMPFR
+
/* Define to 1 if you have the <malloc.h> header file. */
#undef HAVE_MALLOC_H
diff --git a/lib/Makefile.am b/lib/Makefile.am
index 938f7ca1..82459763 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -14,86 +14,89 @@ endif
# Libraries
#
lib_LIBRARIES = libGrid.a
-libGrid_a_SOURCES =\
- Grid_init.cc\
- stencil/Grid_stencil_common.cc\
- qcd/Grid_qcd_dirac.cc\
- qcd/Grid_qcd_wilson_dop.cc\
- algorithms/approx/Zolotarev.cc\
- algorithms/approx/Remez.cc\
+libGrid_a_SOURCES = \
+ Grid_init.cc \
+ stencil/Grid_stencil_common.cc \
+ qcd/Grid_qcd_dirac.cc \
+ qcd/Grid_qcd_wilson_dop.cc \
+ algorithms/approx/Zolotarev.cc \
+ algorithms/approx/Remez.cc \
$(extra_sources)
#
# Include files
#
-nobase_include_HEADERS = algorithms/approx/bigfloat.h\
- algorithms/approx/Chebyshev.h\
- algorithms/approx/Remez.h\
- algorithms/approx/Zolotarev.h\
- algorithms/iterative/ConjugateGradient.h\
- algorithms/iterative/NormalEquations.h\
- algorithms/iterative/SchurRedBlack.h\
- algorithms/LinearOperator.h\
- algorithms/SparseMatrix.h\
- cartesian/Grid_cartesian_base.h\
- cartesian/Grid_cartesian_full.h\
- cartesian/Grid_cartesian_red_black.h\
- communicator/Grid_communicator_base.h\
- cshift/Grid_cshift_common.h\
- cshift/Grid_cshift_mpi.h\
- cshift/Grid_cshift_none.h\
- Grid.h\
- Grid_algorithms.h\
- Grid_aligned_allocator.h\
- Grid_cartesian.h\
- Grid_communicator.h\
- Grid_comparison.h\
- Grid_cshift.h\
- Grid_extract.h\
- Grid_lattice.h\
- Grid_math.h\
- Grid_simd.h\
- Grid_stencil.h\
- Grid_threads.h\
- lattice/Grid_lattice_arith.h\
- lattice/Grid_lattice_base.h\
- lattice/Grid_lattice_comparison.h\
- lattice/Grid_lattice_conformable.h\
- lattice/Grid_lattice_coordinate.h\
- lattice/Grid_lattice_ET.h\
- lattice/Grid_lattice_local.h\
- lattice/Grid_lattice_overload.h\
- lattice/Grid_lattice_peekpoke.h\
- lattice/Grid_lattice_reality.h\
- lattice/Grid_lattice_reduction.h\
- lattice/Grid_lattice_rng.h\
- lattice/Grid_lattice_trace.h\
- lattice/Grid_lattice_transfer.h\
- lattice/Grid_lattice_transpose.h\
- lattice/Grid_lattice_where.h\
- math/Grid_math_arith.h\
- math/Grid_math_arith_add.h\
- math/Grid_math_arith_mac.h\
- math/Grid_math_arith_mul.h\
- math/Grid_math_arith_scalar.h\
- math/Grid_math_arith_sub.h\
- math/Grid_math_inner.h\
- math/Grid_math_outer.h\
- math/Grid_math_peek.h\
- math/Grid_math_poke.h\
- math/Grid_math_reality.h\
- math/Grid_math_tensors.h\
- math/Grid_math_trace.h\
- math/Grid_math_traits.h\
- math/Grid_math_transpose.h\
- parallelIO/GridNerscIO.h\
- qcd/Grid_qcd.h\
- qcd/Grid_qcd_2spinor.h\
- qcd/Grid_qcd_dirac.h\
- qcd/Grid_qcd_wilson_dop.h\
- simd/Grid_vComplexD.h\
- simd/Grid_vComplexF.h\
- simd/Grid_vInteger.h\
- simd/Grid_vRealD.h\
- simd/Grid_vRealF.h
+nobase_include_HEADERS = algorithms/approx/bigfloat.h \
+ algorithms/approx/Chebyshev.h \
+ algorithms/approx/Remez.h \
+ algorithms/approx/Zolotarev.h \
+ algorithms/iterative/ConjugateGradient.h \
+ algorithms/iterative/NormalEquations.h \
+ algorithms/iterative/SchurRedBlack.h \
+ algorithms/LinearOperator.h \
+ algorithms/SparseMatrix.h \
+ cartesian/Grid_cartesian_base.h \
+ cartesian/Grid_cartesian_full.h \
+ cartesian/Grid_cartesian_red_black.h \
+ communicator/Grid_communicator_base.h \
+ cshift/Grid_cshift_common.h \
+ cshift/Grid_cshift_mpi.h \
+ cshift/Grid_cshift_none.h \
+ Grid.h \
+ Grid_algorithms.h \
+ Grid_aligned_allocator.h \
+ Grid_cartesian.h \
+ Grid_communicator.h \
+ Grid_comparison.h \
+ Grid_cshift.h \
+ Grid_extract.h \
+ Grid_lattice.h \
+ Grid_math.h \
+ Grid_simd.h \
+ Grid_stencil.h \
+ Grid_threads.h \
+ lattice/Grid_lattice_arith.h \
+ lattice/Grid_lattice_base.h \
+ lattice/Grid_lattice_comparison.h \
+ lattice/Grid_lattice_conformable.h \
+ lattice/Grid_lattice_coordinate.h \
+ lattice/Grid_lattice_ET.h \
+ lattice/Grid_lattice_local.h \
+ lattice/Grid_lattice_overload.h \
+ lattice/Grid_lattice_peekpoke.h \
+ lattice/Grid_lattice_reality.h \
+ lattice/Grid_lattice_reduction.h \
+ lattice/Grid_lattice_rng.h \
+ lattice/Grid_lattice_trace.h \
+ lattice/Grid_lattice_transfer.h \
+ lattice/Grid_lattice_transpose.h \
+ lattice/Grid_lattice_where.h \
+ math/Grid_math_arith.h \
+ math/Grid_math_arith_add.h \
+ math/Grid_math_arith_mac.h \
+ math/Grid_math_arith_mul.h \
+ math/Grid_math_arith_scalar.h \
+ math/Grid_math_arith_sub.h \
+ math/Grid_math_inner.h \
+ math/Grid_math_outer.h \
+ math/Grid_math_peek.h \
+ math/Grid_math_poke.h \
+ math/Grid_math_reality.h \
+ math/Grid_math_tensors.h \
+ math/Grid_math_trace.h \
+ math/Grid_math_traits.h \
+ math/Grid_math_transpose.h \
+ parallelIO/GridNerscIO.h \
+ qcd/Grid_qcd.h \
+ qcd/Grid_qcd_2spinor.h \
+ qcd/Grid_qcd_dirac.h \
+ qcd/Grid_qcd_wilson_dop.h \
+ simd/Grid_vComplexD.h \
+ simd/Grid_vComplexF.h \
+ simd/Grid_vInteger.h \
+ simd/Grid_vRealD.h \
+ simd/Grid_vRealF.h \
+ simd/Grid_vector_types.h \
+ simd/Grid_sse4.h
+
diff --git a/lib/algorithms/approx/Remez.cc b/lib/algorithms/approx/Remez.cc
old mode 100755
new mode 100644
diff --git a/lib/algorithms/approx/Remez.h b/lib/algorithms/approx/Remez.h
old mode 100755
new mode 100644
diff --git a/lib/algorithms/approx/Zolotarev.cc b/lib/algorithms/approx/Zolotarev.cc
old mode 100755
new mode 100644
diff --git a/lib/algorithms/approx/Zolotarev.h b/lib/algorithms/approx/Zolotarev.h
old mode 100755
new mode 100644
diff --git a/lib/algorithms/approx/bigfloat.h b/lib/algorithms/approx/bigfloat.h
old mode 100755
new mode 100644
diff --git a/lib/simd/.dirstamp b/lib/simd/.dirstamp
new file mode 100644
index 00000000..e69de29b
diff --git a/lib/simd/Grid_sse4.h b/lib/simd/Grid_sse4.h
new file mode 100644
index 00000000..ed4039b7
--- /dev/null
+++ b/lib/simd/Grid_sse4.h
@@ -0,0 +1,194 @@
+//----------------------------------------------------------------------
+/*! @file Grid_sse4.h
+ @brief Optimization libraries
+*/
+// Time-stamp: <2015-05-19 17:06:51 neo>
+//----------------------------------------------------------------------
+
+#include <pmmintrin.h>
+
+namespace Optimization {
+
+ struct Vsplat{
+ //Complex float
+ inline __m128 operator()(float a, float b){
+ return _mm_set_ps(b,a,b,a);
+ }
+ // Real float
+ inline __m128 operator()(float a){
+ return _mm_set_ps(a,a,a,a);
+ }
+ //Complex double
+ inline __m128d operator()(double a, double b){
+ return _mm_set_pd(b,a);
+ }
+ //Real double
+ inline __m128d operator()(double a){
+ return _mm_set_pd(a,a);
+ }
+ //Integer
+ inline __m128i operator()(Integer a){
+ return _mm_set1_epi32(a);
+ }
+ };
+
+ struct Vstore{
+ //Float
+ inline void operator()(__m128 a, float* F){
+ _mm_store_ps(F,a);
+ }
+ //Double
+ inline void operator()(__m128d a, double* D){
+ _mm_store_pd(D,a);
+ }
+ //Integer
+ inline void operator()(__m128i a, Integer* I){
+ _mm_store_si128((__m128i *)I,a);
+ }
+
+ };
+
+
+
+ struct Vset{
+ // Complex float
+ inline __m128 operator()(Grid::ComplexF *a){
+ return _mm_set_ps(a[1].imag(), a[1].real(),a[0].imag(),a[0].real());
+ }
+ // Complex double
+ inline __m128d operator()(Grid::ComplexD *a){
+ return _mm_set_pd(a[0].imag(),a[0].real());
+ }
+ // Real float
+ inline __m128 operator()(float *a){
+ return _mm_set_ps(a[3],a[2],a[1],a[0]);
+ }
+ // Real double
+ inline __m128d operator()(double *a){
+ return _mm_set_pd(a[1],a[0]);
+ }
+ // Integer
+ inline __m128i operator()(Integer *a){
+ return _mm_set_epi32(a[0],a[1],a[2],a[3]);
+ }
+
+
+ };
+
+ struct Reduce{
+ //Complex float
+ inline Grid::ComplexF 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]);
+ }
+ //Complex double
+ inline Grid::ComplexD operator()(__m128d in){
+ printf("Missing complex double implementation -> FIX\n");
+ return Grid::ComplexD(0,0); // FIXME wrong
+ }
+
+
+
+ };
+
+
+ /////////////////////////////////////////////////////
+ // Arithmetic operations
+ /////////////////////////////////////////////////////
+ struct Sum{
+ //Complex/Real float
+ inline __m128 operator()(__m128 a, __m128 b){
+ return _mm_add_ps(a,b);
+ }
+ //Complex/Real double
+ inline __m128d operator()(__m128d a, __m128d b){
+ return _mm_add_pd(a,b);
+ }
+ //Integer
+ inline __m128i operator()(__m128i a, __m128i b){
+ return _mm_add_epi32(a,b);
+ }
+ };
+
+ struct Sub{
+ //Complex/Real float
+ inline __m128 operator()(__m128 a, __m128 b){
+ return _mm_sub_ps(a,b);
+ }
+ //Complex/Real double
+ inline __m128d operator()(__m128d a, __m128d b){
+ return _mm_sub_pd(a,b);
+ }
+ //Integer
+ inline __m128i operator()(__m128i a, __m128i b){
+ return _mm_sub_epi32(a,b);
+ }
+ };
+
+ struct MultComplex{
+ // Complex float
+ inline __m128 operator()(__m128 a, __m128 b){
+ __m128 ymm0,ymm1,ymm2;
+ ymm0 = _mm_shuffle_ps(a,a,_MM_SHUFFLE(2,2,0,0)); // ymm0 <- ar ar,
+ ymm0 = _mm_mul_ps(ymm0,b); // ymm0 <- ar bi, ar br
+ ymm1 = _mm_shuffle_ps(b,b,_MM_SHUFFLE(2,3,0,1)); // ymm1 <- br,bi
+ ymm2 = _mm_shuffle_ps(a,a,_MM_SHUFFLE(3,3,1,1)); // ymm2 <- ai,ai
+ ymm1 = _mm_mul_ps(ymm1,ymm2); // ymm1 <- br ai, ai bi
+ return _mm_addsub_ps(ymm0,ymm1);
+ }
+ // Complex double
+ inline __m128d operator()(__m128d a, __m128d b){
+ __m128d ymm0,ymm1,ymm2;
+ ymm0 = _mm_shuffle_pd(a,a,0x0); // ymm0 <- ar ar,
+ ymm0 = _mm_mul_pd(ymm0,b); // ymm0 <- ar bi, ar br
+ ymm1 = _mm_shuffle_pd(b,b,0x1); // ymm1 <- br,bi b01
+ ymm2 = _mm_shuffle_pd(a,a,0x3); // ymm2 <- ai,ai b11
+ ymm1 = _mm_mul_pd(ymm1,ymm2); // ymm1 <- br ai, ai bi
+ return _mm_addsub_pd(ymm0,ymm1);
+ }
+ };
+
+ struct Mult{
+ // Real float
+ inline __m128 operator()(__m128 a, __m128 b){
+ return _mm_mul_ps(a,b);
+ }
+ // Real double
+ inline __m128d operator()(__m128d a, __m128d b){
+ return _mm_mul_pd(a,b);
+ }
+ // Integer
+ inline __m128i operator()(__m128i a, __m128i b){
+ return _mm_mul_epi32(a,b);
+ }
+
+ };
+
+
+
+
+}
+
+// Here assign types
+namespace Grid {
+ typedef __m128 SIMD_Ftype; // Single precision type
+ typedef __m128d SIMD_Dtype; // Double precision type
+ typedef __m128i SIMD_Itype; // Integer type
+
+
+ // Function names
+ typedef Optimization::Vsplat VsplatSIMD;
+ typedef Optimization::Vstore VstoreSIMD;
+
+ // Arithmetic operations
+ typedef Optimization::Sum SumSIMD;
+ typedef Optimization::Sub SubSIMD;
+ typedef Optimization::Mult MultSIMD;
+ typedef Optimization::MultComplex MultComplexSIMD;
+ typedef Optimization::Vset VsetSIMD;
+
+}
diff --git a/lib/simd/Grid_vector_types.h b/lib/simd/Grid_vector_types.h
index 4576abe3..97958fe4 100644
--- a/lib/simd/Grid_vector_types.h
+++ b/lib/simd/Grid_vector_types.h
@@ -1,6 +1,14 @@
+//---------------------------------------------------------------------------
+/*! @file Grid_vector_types.h
+ @brief Defines templated class Grid_simd to deal with inner vector types
+*/
+// Time-stamp: <2015-05-19 17:20:36 neo>
+//---------------------------------------------------------------------------
#ifndef GRID_VECTOR_TYPES
#define GRID_VECTOR_TYPES
+#include "Grid_sse4.h"
+
namespace Grid {
@@ -13,31 +21,32 @@ namespace Grid {
struct RealPart< std::complex<T> >{
typedef T type;
};
- ////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////
// 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 {};
////////////////////////////////////////////////////////
-
-
// Define the operation templates functors
- template < class SIMD, class Operation >
- SIMD binary(SIMD src_1, SIMD src_2, Operation op){
+ // general forms to allow for vsplat syntax
+ // need explicit declaration of types when used since
+ // clang cannot automatically determine the output type sometimes
+ template < class Out, class Input1, class Input2, class Operation >
+ Out binary(Input1 src_1, Input2 src_2, Operation op){
return op(src_1, src_2);
- }
+ }
- template < class SIMD, class Operation >
- SIMD unary(SIMD src, Operation op){
+ template < class SIMDout, class Input, class Operation >
+ SIMDout unary(Input src, Operation op){
return op(src);
- }
+ }
+
///////////////////////////////////////////////
/*
@brief Grid_simd class for the SIMD vector type operations
-
*/
template < class Scalar_type, class Vector_type >
class Grid_simd {
@@ -69,76 +78,91 @@ namespace Grid {
};
-
///////////////////////////////////////////////
// mac, mult, sub, add, adj
- // Should do an AVX2 version with mac.
///////////////////////////////////////////////
friend inline void mac (Grid_simd * __restrict__ y,const Grid_simd * __restrict__ a,const Grid_simd *__restrict__ x){ *y = (*a)*(*x)+(*y); };
friend inline void mult(Grid_simd * __restrict__ y,const Grid_simd * __restrict__ l,const Grid_simd *__restrict__ r){ *y = (*l) * (*r); }
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); }
+
+ //not for integer types... FIXME
friend inline Grid_simd adj(const Grid_simd &in){ return conjugate(in); }
- //////////////////////////////////
- // Initialise to 1,0,i
- //////////////////////////////////
+ ///////////////////////////////////////////////
+ // Initialise to 1,0,i for the correct types
+ ///////////////////////////////////////////////
// if not complex overload here
- friend inline void vone(Grid_simd &ret) { vsplat(ret,1.0); }
- friend inline void vzero(Grid_simd &ret) { vsplat(ret,0.0); }
-
+ template < class S = Scalar_type,typename std::enable_if < !is_complex < S >::value, int >::type = 0 >
+ friend inline void vone(Grid_simd &ret) { vsplat(ret,1.0); }
+ template < class S = Scalar_type,typename std::enable_if < !is_complex < S >::value, int >::type = 0 >
+ friend inline void vzero(Grid_simd &ret) { vsplat(ret,0.0); }
+
// overload for complex type
template < class S = Scalar_type,typename std::enable_if < is_complex < S >::value, int >::type = 0 >
- friend inline void vone(Grid_simd &ret) { vsplat(ret,1.0,0.0); }
+ friend inline void vone(Grid_simd &ret) { vsplat(ret,1.0,0.0); }
template < class S = Scalar_type,typename std::enable_if < is_complex < S >::value, int >::type = 0 >
- friend inline void vzero(Grid_simd &ret) { vsplat(ret,0.0,0.0); }
-
+ friend inline void vzero(Grid_simd &ret) { vsplat(ret,0.0,0.0); }// use xor?
+
// For integral type
template < class S = Scalar_type,typename std::enable_if < std::is_integral < S >::value, int >::type = 0 >
- friend inline void vone(Grid_simd &ret) { vsplat(ret,1); }
+ friend inline void vone(Grid_simd &ret) { vsplat(ret,1); }
template < class S = Scalar_type,typename std::enable_if < std::is_integral < S >::value, int >::type = 0 >
- friend inline void vzero(Grid_simd &ret) { vsplat(ret,0); }
+ friend inline void vzero(Grid_simd &ret) { vsplat(ret,0); }
+ template < class S = Scalar_type,typename std::enable_if < std::is_integral < S >::value, int >::type = 0 >
+ friend inline void vtrue (Grid_simd &ret){vsplat(ret,0xFFFFFFFF);}
+ template < class S = Scalar_type,typename std::enable_if < std::is_integral < S >::value, int >::type = 0 >
+ friend inline void vfalse(vInteger &ret){vsplat(ret,0);}
-
// do not compile if real or integer, send an error message from the compiler
template < class S = Scalar_type,typename std::enable_if < is_complex < S >::value, int >::type = 0 >
friend inline void vcomplex_i(Grid_simd &ret){ vsplat(ret,0.0,1.0);}
-
-
-
////////////////////////////////////
// Arithmetic operator overloads +,-,*
////////////////////////////////////
friend inline Grid_simd operator + (Grid_simd a, Grid_simd b)
{
- vComplexF ret;
- // FIXME call the binary op
+ 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)
{
- vComplexF ret;
- // FIXME call the binary op
+ Grid_simd ret;
+ ret.v = binary<Vector_type>(a.v, b.v, SubSIMD());
return ret;
};
- friend inline Grid_simd operator * (Grid_simd a, Grid_simd b)
+ // Distinguish between complex types and others
+ template < class S = Scalar_type, typename std::enable_if < is_complex < S >::value, int >::type = 0 >
+ friend inline Grid_simd operator * (Grid_simd a, Grid_simd b)
{
- vComplexF ret;
- // FIXME call the binary op
+ Grid_simd ret;
+ ret.v = binary<Vector_type>(a.v,b.v, MultComplexSIMD());
return ret;
};
-
+
+ // Real/Integer types
+ template < class S = Scalar_type,typename std::enable_if < !is_complex < S >::value, int >::type = 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
////////////////////////////////////////////////////////////////////////
friend inline void vset(Grid_simd &ret, Scalar_type *a){
- // FIXME set
+ ret.v = unary<Vector_type>(a, VsetSIMD());
}
///////////////////////
@@ -147,34 +171,33 @@ namespace Grid {
// overload if complex
template < class S = Scalar_type >
friend inline void vsplat(Grid_simd &ret, typename std::enable_if< is_complex < S >::value, S>::type c){
- Real a= real(c);
- Real b= imag(c);
+ Real a = real(c);
+ Real b = imag(c);
vsplat(ret,a,b);
}
// this only for the complex version
template < class S = Scalar_type, typename std::enable_if < is_complex < S >::value, int >::type = 0 >
friend inline void vsplat(Grid_simd &ret,Real a, Real b){
- // FIXME add operator
+ ret.v = binary<Vector_type>(a, b, VsplatSIMD());
}
- //if real fill with a, if complex fill with a in the real part
+ //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){
- // FIXME add operator
+ ret.v = unary<Vector_type>(a, VsplatSIMD());
}
-
friend inline void vstore(const Grid_simd &ret, Scalar_type *a){
- //FIXME
+ binary<void>(ret.v, (Real*)a, VstoreSIMD());
}
+
friend inline void vprefetch(const Grid_simd &v)
{
_mm_prefetch((const char*)&v.v,_MM_HINT_T0);
}
-
friend inline Scalar_type Reduce(const Grid_simd & in)
{
// FIXME add operator
@@ -221,6 +244,7 @@ namespace Grid {
inline Grid_simd &operator *=(const Grid_simd &r) {
*this = (*this)*r;
return *this;
+ // return (*this)*r; ?
}
inline Grid_simd &operator +=(const Grid_simd &r) {
*this = *this+r;
@@ -233,6 +257,12 @@ namespace Grid {
+ friend inline void permute(Grid_simd &y,Grid_simd b,int perm)
+ {
+ Gpermute<Grid_simd>(y,b,perm);
+ }
+
+ /*
friend inline void permute(Grid_simd &y,Grid_simd b,int perm)
{
Gpermute<Grid_simd>(y,b,perm);
@@ -253,7 +283,7 @@ namespace Grid {
{
Gextract<Grid_simd,Scalar_type>(y,extracted);
}
-
+ */
};// end of Grid_simd class definition
@@ -286,11 +316,11 @@ namespace Grid {
// Define available types (now change names to avoid clashing)
- typedef __m128 SIMD_type;// decided at compilation time
- typedef Grid_simd< float , SIMD_type > MyRealF;
- typedef Grid_simd< double , SIMD_type > MyRealD;
- typedef Grid_simd< std::complex< float > , SIMD_type > MyComplexF;
- typedef Grid_simd< std::complex< double >, SIMD_type > MyComplexD;
+
+ typedef Grid_simd< float , SIMD_Ftype > MyRealF;
+ typedef Grid_simd< double , SIMD_Dtype > MyRealD;
+ typedef Grid_simd< std::complex< float > , SIMD_Ftype > MyComplexF;
+ typedef Grid_simd< std::complex< double >, SIMD_Dtype > MyComplexD;
diff --git a/tests/Grid_main.cc b/tests/Grid_main.cc
index 1aa4fc43..33b1cf4f 100644
--- a/tests/Grid_main.cc
+++ b/tests/Grid_main.cc
@@ -1,5 +1,9 @@
#include "Grid.h"
+//DEBUG
+#include "simd/Grid_vector_types.h"
+
+
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
@@ -151,6 +155,39 @@ int main (int argc, char ** argv)
scMat = sMat*scMat; // LatticeSpinColourMatrix = LatticeSpinMatrix * LatticeSpinColourMatrix
+
+
+#ifdef SSE4
+ ///////// Tests the new class Grid_simd
+ std::complex<double> ctest(3.0,2.0);
+ std::complex<float> ctestf(3.0,2.0);
+ MyComplexF TestMe1(1.0); // fill real part
+ MyComplexD TestMe2(ctest);
+ MyComplexD TestMe3(ctest);// compiler generate conversion of basic types
+ //MyRealF TestMe5(ctest);// Must generate compiler error
+ MyRealD TestMe4(2.0);
+
+ MyComplexF TestMe6(ctestf);
+ MyComplexF TestMe7(ctestf);
+
+ MyComplexD TheSum= TestMe2*TestMe3;
+ MyComplexF TheSumF= TestMe6*TestMe7;
+
+ double dsum[2];
+ _mm_store_pd(dsum, TheSum.v);
+ for (int i =0; i< 2; i++)
+ printf("%f\n", dsum[i]);
+
+ float fsum[4];
+ _mm_store_ps(fsum, TheSumF.v);
+ for (int i =0; i< 4; i++)
+ printf("%f\n", fsum[i]);
+
+ vstore(TheSum, &ctest);
+ std::cout << ctest<< std::endl;
+#endif
+ ///////////////////////
+
// Non-lattice (const objects) * Lattice
ColourMatrix cm;
SpinColourMatrix scm;