Merge branch 'develop' of https://github.com/paboyle/Grid into feature/staggering

.gitignore

@@ -47,7 +47,8 @@ Config.h.in
 config.log
 config.status
 .deps
-*.inc
+Make.inc
+eigen.inc
 
 # http://www.gnu.org/software/autoconf #
 ########################################

Makefile.am

@@ -1,10 +1,12 @@
 # additional include paths necessary to compile the C++ library
 SUBDIRS = lib benchmarks tests
 
-.PHONY: tests
+include $(top_srcdir)/doxygen.inc
 
 tests: all
 	$(MAKE) -C tests tests
 
+.PHONY: tests doxygen-run doxygen-doc $(DX_PS_GOAL) $(DX_PDF_GOAL)
+
 AM_CXXFLAGS += -I$(top_builddir)/include
 ACLOCAL_AMFLAGS = -I m4

README.md

@@ -116,13 +116,15 @@ If you want to build all the tests at once just use `make tests`.
 - `--with-fftw=<path>`: look for FFTW in the UNIX prefix `<path>`
 - `--enable-lapack[=<path>]`: enable LAPACK support in Lanczos eigensolver. A UNIX prefix containing the library can be specified (optional).
 - `--enable-mkl[=<path>]`: use Intel MKL for FFT (and LAPACK if enabled) routines. A UNIX prefix containing the library can be specified (optional).
-- `--enable-numa`: ???
+- `--enable-numa`: enable NUMA first touch optimisation
 - `--enable-simd=<code>`: setup Grid for the SIMD target `<code>` (default: `GEN`). A list of possible SIMD targets is detailed in a section below.
+- `--enable-gen-simd-width=<size>`: select the size (in bytes) of the generic SIMD vector type (default: 32 bytes).
 - `--enable-precision={single|double}`: set the default precision (default: `double`).
 - `--enable-precision=<comm>`: Use `<comm>` for message passing (default: `none`). A list of possible SIMD targets is detailed in a section below.
 - `--enable-rng={ranlux48|mt19937}`: choose the RNG (default: `ranlux48 `).
 - `--disable-timers`: disable system dependent high-resolution timers.
 - `--enable-chroma`: enable Chroma regression tests.
+- `--enable-doxygen-doc`: enable the Doxygen documentation generation (build with `make doxygen-doc`)
 
 ### Possible communication interfaces
 
@@ -165,6 +167,7 @@ Alternatively, some CPU codenames can be directly used:
 - We currently support AVX512 only for the Intel compiler. Support for GCC and clang will appear in future versions of Grid when the AVX512 support within GCC and clang will be more advanced.
 - For BG/Q only [bgclang](http://trac.alcf.anl.gov/projects/llvm-bgq) is supported. We do not presently plan to support more compilers for this platform.
 - BG/Q performances are currently rather poor. This is being investigated for future versions.
+- The vector size for the `GEN` target can be specified with the `configure` script option `--enable-gen-simd-width`.
 
 ### Build setup for Intel Knights Landing platform
 

benchmarks/Makefile.am

@@ -1 +1,11 @@
 include Make.inc
+
+simple: simple_su3_test.o simple_su3_expr.o simple_simd_test.o
+
+EXTRA_LIBRARIES = libsimple_su3_test.a libsimple_su3_expr.a libsimple_simd_test.a
+
+libsimple_su3_test_a_SOURCES = simple_su3_test.cc
+
+libsimple_su3_expr_a_SOURCES = simple_su3_expr.cc
+
+libsimple_simd_test_a_SOURCES = simple_simd_test.cc

benchmarks/simple_simd_test.cc

@@ -0,0 +1,11 @@
+#include <Grid/Grid.h>
+
+Grid::vRealD add(const Grid::vRealD &x, const Grid::vRealD &y)
+{
+	return x+y;
+}
+
+Grid::vRealD sub(const Grid::vRealD &x, const Grid::vRealD &y)
+{
+	return x-y;
+}

benchmarks/simple_su3_expr.cc

@@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 using namespace std;
 using namespace Grid;

benchmarks/simple_su3_test.cc

@@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include <Grid.h>
+#include <Grid/Grid.h>
 
 using namespace std;
 using namespace Grid;

configure.ac

@@ -149,8 +149,14 @@ CXXFLAGS=$CXXFLAGS_CPY
 LDFLAGS=$LDFLAGS_CPY
 
 ############### SIMD instruction selection
-AC_ARG_ENABLE([simd],[AC_HELP_STRING([--enable-simd=<code>],
-	            [select SIMD target (cf. README.md)])], [ac_SIMD=${enable_simd}], [ac_SIMD=GEN])
+AC_ARG_ENABLE([simd],[AC_HELP_STRING([--enable-simd=code],
+              [select SIMD target (cf. README.md)])], [ac_SIMD=${enable_simd}], [ac_SIMD=GEN])
+
+AC_ARG_ENABLE([gen-simd-width],
+            [AS_HELP_STRING([--enable-gen-simd-width=size],
+            [size (in bytes) of the generic SIMD vectors (default: 32)])],
+            [ac_gen_simd_width=$enable_gen_simd_width],
+            [ac_gen_simd_width=32])
 
 case ${ax_cv_cxx_compiler_vendor} in
   clang|gnu)
@@ -180,7 +186,10 @@ case ${ax_cv_cxx_compiler_vendor} in
         AC_DEFINE([AVX512],[1],[AVX512 intrinsics])
         SIMD_FLAGS='-march=knl';;
       GEN)
-        AC_DEFINE([GENERIC_VEC],[1],[generic vector code])
+        AC_DEFINE([GEN],[1],[generic vector code])
+        AC_DEFINE_UNQUOTED([GEN_SIMD_WIDTH],[$ac_gen_simd_width],
+                           [generic SIMD vector width (in bytes)])
+        SIMD_GEN_WIDTH_MSG=" (width= $ac_gen_simd_width)"
         SIMD_FLAGS='';;
       QPX|BGQ)
         AC_DEFINE([QPX],[1],[QPX intrinsics for BG/Q])
@@ -212,7 +221,10 @@ case ${ax_cv_cxx_compiler_vendor} in
         AC_DEFINE([AVX512],[1],[AVX512 intrinsics for Knights Landing])
         SIMD_FLAGS='-xmic-avx512';;
       GEN)
-        AC_DEFINE([GENERIC_VEC],[1],[generic vector code])
+        AC_DEFINE([GEN],[1],[generic vector code])
+        AC_DEFINE_UNQUOTED([GEN_SIMD_WIDTH],[$ac_gen_simd_width],
+                           [generic SIMD vector width (in bytes)])
+        SIMD_GEN_WIDTH_MSG=" (width= $ac_gen_simd_width)"
         SIMD_FLAGS='';;
       *)
         AC_MSG_ERROR(["SIMD option ${ac_SIMD} not supported by the Intel compiler"]);;
@@ -342,12 +354,17 @@ esac
 AM_CONDITIONAL(BUILD_CHROMA_REGRESSION,[ test "X${ac_CHROMA}X" == "XyesX" ])
 
 ############### Doxygen
-AC_PROG_DOXYGEN
-
-if test -n "$DOXYGEN"
-then
-AC_CONFIG_FILES([docs/doxy.cfg])
-fi
+DX_DOXYGEN_FEATURE([OFF])
+DX_DOT_FEATURE([OFF])
+DX_HTML_FEATURE([ON])
+DX_CHM_FEATURE([OFF])
+DX_CHI_FEATURE([OFF])
+DX_MAN_FEATURE([OFF])
+DX_RTF_FEATURE([OFF])
+DX_XML_FEATURE([OFF])
+DX_PDF_FEATURE([OFF])
+DX_PS_FEATURE([OFF])
+DX_INIT_DOXYGEN([$PACKAGE_NAME], [doxygen.cfg])
 
 ############### Ouput
 cwd=`pwd -P`; cd ${srcdir}; abs_srcdir=`pwd -P`; cd ${cwd}
@@ -382,7 +399,7 @@ os (target)                 : $target_os
 compiler vendor             : ${ax_cv_cxx_compiler_vendor}
 compiler version            : ${ax_cv_gxx_version}
 ----- BUILD OPTIONS -----------------------------------
-SIMD                        : ${ac_SIMD}
+SIMD                        : ${ac_SIMD}${SIMD_GEN_WIDTH_MSG}
 Threading                   : ${ac_openmp} 
 Communications type         : ${comms_type}
 Default precision           : ${ac_PRECISION}
@@ -390,8 +407,7 @@ RNG choice                  : ${ac_RNG}
 GMP                         : `if test "x$have_gmp" = xtrue; then echo yes; else echo no; fi`
 LAPACK                      : ${ac_LAPACK}
 FFTW                        : `if test "x$have_fftw" = xtrue; then echo yes; else echo no; fi`
-build DOXYGEN documentation : `if test "x$enable_doc" = xyes; then echo yes; else echo no; fi`
-graphs and diagrams         : `if test "x$enable_dot" = xyes; then echo yes; else echo no; fi`
+build DOXYGEN documentation : `if test "$DX_FLAG_doc" = '1'; then echo yes; else echo no; fi`
 ----- BUILD FLAGS -------------------------------------
 CXXFLAGS:
 `echo ${AM_CXXFLAGS} ${CXXFLAGS} | tr ' ' '\n' | sed 's/^-/    -/g'`

doxygen.inc

@@ -0,0 +1,184 @@
+# Copyright (C) 2004 Oren Ben-Kiki
+# This file is distributed under the same terms as the Automake macro files.
+
+# Generate automatic documentation using Doxygen. Goals and variables values
+# are controlled by the various DX_COND_??? conditionals set by autoconf.
+#
+# The provided goals are:
+# doxygen-doc: Generate all doxygen documentation.
+# doxygen-run: Run doxygen, which will generate some of the documentation
+#              (HTML, CHM, CHI, MAN, RTF, XML) but will not do the post
+#              processing required for the rest of it (PS, PDF, and some MAN).
+# doxygen-man: Rename some doxygen generated man pages.
+# doxygen-ps: Generate doxygen PostScript documentation.
+# doxygen-pdf: Generate doxygen PDF documentation.
+#
+# Note that by default these are not integrated into the automake goals. If
+# doxygen is used to generate man pages, you can achieve this integration by
+# setting man3_MANS to the list of man pages generated and then adding the
+# dependency:
+#
+#   $(man3_MANS): doxygen-doc
+#
+# This will cause make to run doxygen and generate all the documentation.
+#
+# The following variable is intended for use in Makefile.am:
+#
+# DX_CLEANFILES = everything to clean.
+#
+# This is usually added to MOSTLYCLEANFILES.
+
+## --------------------------------- ##
+## Format-independent Doxygen rules. ##
+## --------------------------------- ##
+
+if DX_COND_doc
+
+## ------------------------------- ##
+## Rules specific for HTML output. ##
+## ------------------------------- ##
+
+if DX_COND_html
+
+DX_CLEAN_HTML = @DX_DOCDIR@/html
+
+endif DX_COND_html
+
+## ------------------------------ ##
+## Rules specific for CHM output. ##
+## ------------------------------ ##
+
+if DX_COND_chm
+
+DX_CLEAN_CHM = @DX_DOCDIR@/chm
+
+if DX_COND_chi
+
+DX_CLEAN_CHI = @DX_DOCDIR@/@PACKAGE@.chi
+
+endif DX_COND_chi
+
+endif DX_COND_chm
+
+## ------------------------------ ##
+## Rules specific for MAN output. ##
+## ------------------------------ ##
+
+if DX_COND_man
+
+DX_CLEAN_MAN = @DX_DOCDIR@/man
+
+endif DX_COND_man
+
+## ------------------------------ ##
+## Rules specific for RTF output. ##
+## ------------------------------ ##
+
+if DX_COND_rtf
+
+DX_CLEAN_RTF = @DX_DOCDIR@/rtf
+
+endif DX_COND_rtf
+
+## ------------------------------ ##
+## Rules specific for XML output. ##
+## ------------------------------ ##
+
+if DX_COND_xml
+
+DX_CLEAN_XML = @DX_DOCDIR@/xml
+
+endif DX_COND_xml
+
+## ----------------------------- ##
+## Rules specific for PS output. ##
+## ----------------------------- ##
+
+if DX_COND_ps
+
+DX_CLEAN_PS = @DX_DOCDIR@/@PACKAGE@.ps
+
+DX_PS_GOAL = doxygen-ps
+
+doxygen-ps: @DX_DOCDIR@/@PACKAGE@.ps
+
+@DX_DOCDIR@/@PACKAGE@.ps: @DX_DOCDIR@/@PACKAGE@.tag
+	cd @DX_DOCDIR@/latex; \
+	rm -f *.aux *.toc *.idx *.ind *.ilg *.log *.out; \
+	$(DX_LATEX) refman.tex; \
+	$(MAKEINDEX_PATH) refman.idx; \
+	$(DX_LATEX) refman.tex; \
+	countdown=5; \
+	while $(DX_EGREP) 'Rerun (LaTeX|to get cross-references right)' \
+	                  refman.log > /dev/null 2>&1 \
+	   && test $$countdown -gt 0; do \
+	    $(DX_LATEX) refman.tex; \
+	    countdown=`expr $$countdown - 1`; \
+	done; \
+	$(DX_DVIPS) -o ../@PACKAGE@.ps refman.dvi
+
+endif DX_COND_ps
+
+## ------------------------------ ##
+## Rules specific for PDF output. ##
+## ------------------------------ ##
+
+if DX_COND_pdf
+
+DX_CLEAN_PDF = @DX_DOCDIR@/@PACKAGE@.pdf
+
+DX_PDF_GOAL = doxygen-pdf
+
+doxygen-pdf: @DX_DOCDIR@/@PACKAGE@.pdf
+
+@DX_DOCDIR@/@PACKAGE@.pdf: @DX_DOCDIR@/@PACKAGE@.tag
+	cd @DX_DOCDIR@/latex; \
+	rm -f *.aux *.toc *.idx *.ind *.ilg *.log *.out; \
+	$(DX_PDFLATEX) refman.tex; \
+	$(DX_MAKEINDEX) refman.idx; \
+	$(DX_PDFLATEX) refman.tex; \
+	countdown=5; \
+	while $(DX_EGREP) 'Rerun (LaTeX|to get cross-references right)' \
+	                  refman.log > /dev/null 2>&1 \
+	   && test $$countdown -gt 0; do \
+	    $(DX_PDFLATEX) refman.tex; \
+	    countdown=`expr $$countdown - 1`; \
+	done; \
+	mv refman.pdf ../@PACKAGE@.pdf
+
+endif DX_COND_pdf
+
+## ------------------------------------------------- ##
+## Rules specific for LaTeX (shared for PS and PDF). ##
+## ------------------------------------------------- ##
+
+if DX_COND_latex
+
+DX_CLEAN_LATEX = @DX_DOCDIR@/latex
+
+endif DX_COND_latex
+
+.INTERMEDIATE: doxygen-run $(DX_PS_GOAL) $(DX_PDF_GOAL)
+
+doxygen-run: @DX_DOCDIR@/@PACKAGE@.tag
+
+doxygen-doc: doxygen-run $(DX_PS_GOAL) $(DX_PDF_GOAL)
+
+@DX_DOCDIR@/@PACKAGE@.tag: $(DX_CONFIG) $(pkginclude_HEADERS)
+	rm -rf @DX_DOCDIR@
+	$(DX_ENV) $(DX_DOXYGEN) $(srcdir)/$(DX_CONFIG)
+
+DX_CLEANFILES = \
+    @DX_DOCDIR@/@PACKAGE@.tag \
+    -r \
+    $(DX_CLEAN_HTML) \
+    $(DX_CLEAN_CHM) \
+    $(DX_CLEAN_CHI) \
+    $(DX_CLEAN_MAN) \
+    $(DX_CLEAN_RTF) \
+    $(DX_CLEAN_XML) \
+    $(DX_CLEAN_PS) \
+    $(DX_CLEAN_PDF) \
+    $(DX_CLEAN_LATEX)
+
+endif DX_COND_doc

lib/algorithms/iterative/ImplicitlyRestartedLanczos.h

@@ -1080,10 +1080,10 @@ say con = 2
 **/
 
 template<class T>
-static void Lock(DenseMatrix<T> &H, 	///Hess mtx	
-		 DenseMatrix<T> &Q, 	///Lock Transform
-		 T val, 		///value to be locked
-		 int con, 	///number already locked
+static void Lock(DenseMatrix<T> &H, 	// Hess mtx	
+		 DenseMatrix<T> &Q, 	// Lock Transform
+		 T val, 		// value to be locked
+		 int con, 	// number already locked
 		 RealD small,
 		 int dfg,
 		 bool herm)

lib/simd/Grid_avx512.h

@@ -27,15 +27,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-//----------------------------------------------------------------------
-/*! @file Grid_knc.h
-  @brief Optimization libraries for AVX512 instructions set for KNC
-
-  Using intrinsics
-*/
-// Time-stamp: <2015-06-09 14:27:28 neo>
-//----------------------------------------------------------------------
-
 #include <immintrin.h>
 
 

lib/simd/Grid_generic.h

@@ -6,8 +6,7 @@
 
     Copyright (C) 2015
 
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: neo <cossu@post.kek.jp>
+Author: Antonin Portelli <antonin.portelli@me.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
@@ -27,133 +26,352 @@ Author: neo <cossu@post.kek.jp>
     *************************************************************************************/
     /*  END LEGAL */
 
+static_assert(GEN_SIMD_WIDTH % 16u == 0, "SIMD vector size is not an integer multiple of 16 bytes");
+
+//#define VECTOR_LOOPS
+
+// playing with compiler pragmas
+#ifdef VECTOR_LOOPS
+#ifdef __clang__
+#define VECTOR_FOR(i, w, inc)\
+_Pragma("clang loop unroll(full) vectorize(enable) interleave(enable) vectorize_width(w)")\
+for (unsigned int i = 0; i < w; i += inc)
+#elif defined __INTEL_COMPILER
+#define VECTOR_FOR(i, w, inc)\
+_Pragma("simd vectorlength(w*8)")\
+for (unsigned int i = 0; i < w; i += inc)
+#else
+#define VECTOR_FOR(i, w, inc)\
+for (unsigned int i = 0; i < w; i += inc)
+#endif
+#else
+#define VECTOR_FOR(i, w, inc)\
+for (unsigned int i = 0; i < w; i += inc)
+#endif
+
 namespace Grid {
 namespace Optimization {
 
-  template<class vtype>
-  union uconv {
-    float f;
-    vtype v;
+  // type traits giving the number of elements for each vector type
+  template <typename T> struct W;
+  template <> struct W<double> {
+    constexpr static unsigned int c = GEN_SIMD_WIDTH/16u;
+    constexpr static unsigned int r = GEN_SIMD_WIDTH/8u;
   };
-
-  union u128f {
-    float v;
-    float f[4];
-  };
-  union u128d {
-    double v;
-    double f[2];
+  template <> struct W<float> {
+    constexpr static unsigned int c = GEN_SIMD_WIDTH/8u;
+    constexpr static unsigned int r = GEN_SIMD_WIDTH/4u;
   };
   
+  // SIMD vector types
+  template <typename T>
+  struct vec {
+    alignas(GEN_SIMD_WIDTH) T v[W<T>::r];
+  };
+  
+  typedef vec<float>   vecf;
+  typedef vec<double>  vecd;
+  
   struct Vsplat{
-    //Complex float
-    inline u128f operator()(float a, float b){
-      u128f out; 
-      out.f[0] = a;
-      out.f[1] = b;
-      out.f[2] = a;
-      out.f[3] = b;
+    // Complex
+    template <typename T>
+    inline vec<T> operator()(T a, T b){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::r, 2)
+      {
+        out.v[i]   = a;
+        out.v[i+1] = b;
+      }
+
       return out;
     }
-    // Real float
-    inline u128f operator()(float a){
-      u128f out; 
-      out.f[0] = a;
-      out.f[1] = a;
-      out.f[2] = a;
-      out.f[3] = a;
+    
+    // Real
+    template <typename T>
+    inline vec<T> operator()(T a){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::r, 1)
+      {
+        out.v[i] = a;
+      }
+      
       return out;
     }
-    //Complex double
-    inline u128d operator()(double a, double b){
-      u128d out; 
-      out.f[0] = a;
-      out.f[1] = b;
-      return out;
-    }
-    //Real double
-    inline u128d operator()(double a){
-      u128d out; 
-      out.f[0] = a;
-      out.f[1] = a;
-      return out;
-    }
-    //Integer
+    
+    // Integer
     inline int operator()(Integer a){
       return a;
     }
   };
 
   struct Vstore{
-    //Float 
-    inline void operator()(u128f a, float* F){
-      memcpy(F,a.f,4*sizeof(float));
-    }
-    //Double
-    inline void operator()(u128d a, double* D){
-      memcpy(D,a.f,2*sizeof(double));
+    // Real
+    template <typename T>
+    inline void operator()(vec<T> a, T *D){
+      *((vec<T> *)D) = a;
     }
     //Integer
-    inline void operator()(int a, Integer* I){
-      I[0] = a;
+    inline void operator()(int a, Integer *I){
+      *I = a;
     }
 
   };
 
   struct Vstream{
-    //Float
-    inline void operator()(float * a, u128f b){
-      memcpy(a,b.f,4*sizeof(float));
+    // Real
+    template <typename T>
+    inline void operator()(T * a, vec<T> b){
+      *((vec<T> *)a) = b;
     }
-    //Double
-    inline void operator()(double * a, u128d b){
-      memcpy(a,b.f,2*sizeof(double));
-    }
-
-
   };
 
   struct Vset{
-    // Complex float 
-    inline u128f operator()(Grid::ComplexF *a){
-      u128f out; 
-      out.f[0] = a[0].real();
-      out.f[1] = a[0].imag();
-      out.f[2] = a[1].real();
-      out.f[3] = a[1].imag();
+    // Complex
+    template <typename T>
+    inline vec<T> operator()(std::complex<T> *a){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::c, 1)
+      {
+        out.v[2*i]   = a[i].real();
+        out.v[2*i+1] = a[i].imag();
+      }
+      
       return out;
     }
-    // Complex double 
-    inline u128d operator()(Grid::ComplexD *a){
-      u128d out; 
-      out.f[0] = a[0].real();
-      out.f[1] = a[0].imag();
-      return out;
-    }
-    // Real float 
-    inline u128f operator()(float *a){
-      u128f out; 
-      out.f[0] = a[0];
-      out.f[1] = a[1];
-      out.f[2] = a[2];
-      out.f[3] = a[3];
-      return out;
-    }
-    // Real double
-    inline u128d operator()(double *a){
-      u128d out; 
-      out.f[0] = a[0];
-      out.f[1] = a[1];
+    
+    // Real
+    template <typename T>
+    inline vec<T> operator()(T *a){
+      vec<T> out;
+      
+      out = *((vec<T> *)a);
+      
       return out;
     }
+
     // Integer
     inline int operator()(Integer *a){
-      return a[0];
+      return *a;
     }
-
-
   };
 
+  /////////////////////////////////////////////////////
+  // Arithmetic operations
+  /////////////////////////////////////////////////////
+  struct Sum{
+    // Complex/Real
+    template <typename T>
+    inline vec<T> operator()(vec<T> a, vec<T> b){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::r, 1)
+      {
+        out.v[i] = a.v[i] + b.v[i];
+      }
+      
+      return out;
+    }
+    
+    //I nteger
+    inline int operator()(int a, int b){
+      return a + b;
+    }
+  };
+
+  struct Sub{
+    // Complex/Real
+    template <typename T>
+    inline vec<T> operator()(vec<T> a, vec<T> b){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::r, 1)
+      {
+        out.v[i] = a.v[i] - b.v[i];
+      }
+      
+      return out;
+    }
+    
+    //Integer
+    inline int operator()(int a, int b){
+      return a-b;
+    }
+  };
+
+  struct Mult{
+    // Real
+    template <typename T>
+    inline vec<T> operator()(vec<T> a, vec<T> b){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::r, 1)
+      {
+        out.v[i] = a.v[i]*b.v[i];
+      }
+      
+      return out;
+    }
+    
+    // Integer
+    inline int operator()(int a, int b){
+      return a*b;
+    }
+  };
+  
+  #define cmul(a, b, c, i)\
+  c[i]   = a[i]*b[i]   - a[i+1]*b[i+1];\
+  c[i+1] = a[i]*b[i+1] + a[i+1]*b[i];
+  
+  struct MultComplex{
+    // Complex
+    template <typename T>
+    inline vec<T> operator()(vec<T> a, vec<T> b){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::c, 1)
+      {
+        cmul(a.v, b.v, out.v, 2*i);
+      }      
+      
+      return out;
+    }
+  };
+  
+  #undef cmul
+
+  struct Div{
+    // Real
+    template <typename T>
+    inline vec<T> operator()(vec<T> a, vec<T> b){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::r, 1)
+      {
+        out.v[i] = a.v[i]/b.v[i];
+      }
+      
+      return out;
+    }
+  };
+  
+  #define conj(a, b, i)\
+  b[i]   = a[i];\
+  b[i+1] = -a[i+1];
+  
+  struct Conj{
+    // Complex
+    template <typename T>
+    inline vec<T> operator()(vec<T> a){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::c, 1)
+      {
+        conj(a.v, out.v, 2*i);
+      }
+      
+      return out;
+    }
+  };
+  
+  #undef conj
+
+  #define timesmi(a, b, i)\
+  b[i]   = a[i+1];\
+  b[i+1] = -a[i];
+  
+  struct TimesMinusI{
+    // Complex
+    template <typename T>
+    inline vec<T> operator()(vec<T> a, vec<T> b){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::c, 1)
+      {
+        timesmi(a.v, out.v, 2*i);
+      }
+      
+      return out;
+    }
+  };
+
+  #undef timesmi
+  
+  #define timesi(a, b, i)\
+  b[i]   = -a[i+1];\
+  b[i+1] = a[i];
+  
+  struct TimesI{
+    // Complex
+    template <typename T>
+    inline vec<T> operator()(vec<T> a, vec<T> b){
+      vec<T> out;
+      
+      VECTOR_FOR(i, W<T>::c, 1)
+      {
+        timesi(a.v, out.v, 2*i);
+      }
+      
+      return out;
+    }
+  };
+  
+  #undef timesi
+
+  //////////////////////////////////////////////
+  // Some Template specialization
+  #define perm(a, b, n, w)\
+  unsigned int _mask = w >> (n + 1);\
+  VECTOR_FOR(i, w, 1)\
+  {\
+    b[i] = a[i^_mask];\
+  }
+  
+  #define DECL_PERMUTE_N(n)\
+  template <typename T>\
+  static inline vec<T> Permute##n(vec<T> in) {\
+    vec<T> out;\
+    perm(in.v, out.v, n, W<T>::r);\
+    return out;\
+  }
+  
+  struct Permute{
+    DECL_PERMUTE_N(0);
+    DECL_PERMUTE_N(1);
+    DECL_PERMUTE_N(2);
+    DECL_PERMUTE_N(3);
+  };
+  
+  #undef perm
+  #undef DECL_PERMUTE_N
+  
+  #define rot(a, b, n, w)\
+  VECTOR_FOR(i, w, 1)\
+  {\
+    b[i] = a[(i + n)%w];\
+  }
+  
+  struct Rotate{
+    template <typename T>
+    static inline vec<T> rotate(vec<T> in, int n){
+      vec<T> out;
+      
+      rot(in.v, out.v, n, W<T>::r);
+      
+      return out;
+    }
+  };
+
+  #undef rot
+  
+  #define acc(v, a, off, step, n)\
+  for (unsigned int i = off; i < n; i += step)\
+  {\
+    a += v[i];\
+  }
+  
   template <typename Out_type, typename In_type>
   struct Reduce{
     //Need templated class to overload output type
@@ -164,316 +382,67 @@ namespace Optimization {
       return 0;
     }
   };
-
-  /////////////////////////////////////////////////////
-  // Arithmetic operations
-  /////////////////////////////////////////////////////
-  struct Sum{
-    //Complex/Real float
-    inline u128f operator()(u128f a, u128f b){
-      u128f out;
-      out.f[0] = a.f[0] + b.f[0];
-      out.f[1] = a.f[1] + b.f[1];
-      out.f[2] = a.f[2] + b.f[2];
-      out.f[3] = a.f[3] + b.f[3];
-      return out;
-    }
-    //Complex/Real double
-    inline u128d operator()(u128d a, u128d b){
-      u128d out;
-      out.f[0] = a.f[0] + b.f[0];
-      out.f[1] = a.f[1] + b.f[1];
-      return out;
-    }
-    //Integer
-    inline int operator()(int a, int b){
-      return a + b;
-    }
-  };
-
-  struct Sub{
-    //Complex/Real float
-    inline u128f operator()(u128f a, u128f b){
-      u128f out;
-      out.f[0] = a.f[0] - b.f[0];
-      out.f[1] = a.f[1] - b.f[1];
-      out.f[2] = a.f[2] - b.f[2];
-      out.f[3] = a.f[3] - b.f[3];
-      return out;
-    }
-    //Complex/Real double
-    inline u128d operator()(u128d a, u128d b){
-      u128d out;
-      out.f[0] = a.f[0] - b.f[0];
-      out.f[1] = a.f[1] - b.f[1];
-      return out;
-    }
-    //Integer
-    inline int operator()(int a, int b){
-      return a-b;
-    }
-  };
-
-  struct MultComplex{
-    // Complex float
-    inline u128f operator()(u128f a, u128f b){
-      u128f out;
-      out.f[0] = a.f[0]*b.f[0] - a.f[1]*b.f[1];
-      out.f[1] = a.f[0]*b.f[1] + a.f[1]*b.f[0];
-      out.f[2] = a.f[2]*b.f[2] - a.f[3]*b.f[3];
-      out.f[3] = a.f[2]*b.f[3] + a.f[3]*b.f[2];
-      return out;
-    }
-    // Complex double
-    inline u128d operator()(u128d a, u128d b){
-      u128d out;
-      out.f[0] = a.f[0]*b.f[0] - a.f[1]*b.f[1];
-      out.f[1] = a.f[0]*b.f[1] + a.f[1]*b.f[0];
-      return out;
-    }
-  };
-
-  struct Mult{
-    //CK: Appear unneeded
-    // inline float  mac(float a, float b,double c){
-    //   return 0;
-    // }
-    // inline double mac(double a, double b,double c){
-    //   return 0;
-    // }
-
-    // Real float
-    inline u128f operator()(u128f a, u128f b){
-      u128f out;
-      out.f[0] = a.f[0]*b.f[0];
-      out.f[1] = a.f[1]*b.f[1];
-      out.f[2] = a.f[2]*b.f[2];
-      out.f[3] = a.f[3]*b.f[3];
-      return out;
-    }
-    // Real double
-    inline u128d operator()(u128d a, u128d b){
-      u128d out;
-      out.f[0] = a.f[0]*b.f[0];
-      out.f[1] = a.f[1]*b.f[1];
-      return out;
-    }
-    // Integer
-    inline int operator()(int a, int b){
-      return a*b;
-    }
-  };
-
-  struct Conj{
-    // Complex single
-    inline u128f operator()(u128f in){
-      u128f out;
-      out.f[0] = in.f[0];
-      out.f[1] = -in.f[1];
-      out.f[2] = in.f[2];
-      out.f[3] = -in.f[3];
-      return out;
-    }
-    // Complex double
-    inline u128d operator()(u128d in){
-      u128d out;
-      out.f[0] = in.f[0];
-      out.f[1] = -in.f[1];
-      return out;
-    }
-    // do not define for integer input
-  };
-
-  struct TimesMinusI{
-    //Complex single
-    inline u128f operator()(u128f in, u128f ret){ //note ret is ignored
-      u128f out;
-      out.f[0] = in.f[1];
-      out.f[1] = -in.f[0];
-      out.f[2] = in.f[3];
-      out.f[3] = -in.f[2];
-      return out;
-    }
-    //Complex double
-    inline u128d operator()(u128d in, u128d ret){
-      u128d out;
-      out.f[0] = in.f[1];
-      out.f[1] = -in.f[0];
-      return out;
-    }
-  };
-
-  struct TimesI{
-    //Complex single
-    inline u128f operator()(u128f in, u128f ret){ //note ret is ignored
-      u128f out;
-      out.f[0] = -in.f[1];
-      out.f[1] = in.f[0];
-      out.f[2] = -in.f[3];
-      out.f[3] = in.f[2];
-      return out;
-    }
-    //Complex double
-    inline u128d operator()(u128d in, u128d ret){
-      u128d out;
-      out.f[0] = -in.f[1];
-      out.f[1] = in.f[0];
-      return out;
-    }
-  };
-
-  //////////////////////////////////////////////
-  // Some Template specialization
-  struct Permute{
-    //We just have to mirror the permutes of Grid_sse4.h
-    static inline u128f Permute0(u128f in){ //AB CD -> CD AB
-      u128f out;
-      out.f[0] = in.f[2];
-      out.f[1] = in.f[3];
-      out.f[2] = in.f[0];
-      out.f[3] = in.f[1];
-      return out;
-    };
-    static inline u128f Permute1(u128f in){ //AB CD -> BA DC
-      u128f out;
-      out.f[0] = in.f[1];
-      out.f[1] = in.f[0];
-      out.f[2] = in.f[3];
-      out.f[3] = in.f[2];
-      return out;
-    };
-    static inline u128f Permute2(u128f in){
-      return in;
-    };
-    static inline u128f Permute3(u128f in){
-      return in;
-    };
-
-    static inline u128d Permute0(u128d in){ //AB -> BA
-      u128d out;
-      out.f[0] = in.f[1];
-      out.f[1] = in.f[0];
-      return out;      
-    };
-    static inline u128d Permute1(u128d in){
-      return in;
-    };
-    static inline u128d Permute2(u128d in){
-      return in;
-    };
-    static inline u128d Permute3(u128d in){
-      return in;
-    };
-
-  };
   
-  template < typename vtype > 
-    void permute(vtype &a, vtype b, int perm) {
-   };
-    
-  struct Rotate{
-
-    static inline u128f rotate(u128f in,int n){
-      u128f out;
-      switch(n){
-      case 0:
-        out.f[0] = in.f[0];
-        out.f[1] = in.f[1];
-        out.f[2] = in.f[2];
-        out.f[3] = in.f[3];
-        break;
-      case 1:
-        out.f[0] = in.f[1];
-        out.f[1] = in.f[2];
-        out.f[2] = in.f[3];
-        out.f[3] = in.f[0];
-        break;
-      case 2:
-        out.f[0] = in.f[2];
-        out.f[1] = in.f[3];
-        out.f[2] = in.f[0];
-        out.f[3] = in.f[1];
-        break;
-      case 3:
-        out.f[0] = in.f[3];
-        out.f[1] = in.f[0];
-        out.f[2] = in.f[1];
-        out.f[3] = in.f[2];
-        break;
-      default: assert(0);
-      }
-      return out;
-    }
-    static inline u128d rotate(u128d in,int n){
-      u128d out;
-      switch(n){
-      case 0:
-        out.f[0] = in.f[0];
-        out.f[1] = in.f[1];
-        break;
-      case 1:
-        out.f[0] = in.f[1];
-        out.f[1] = in.f[0];
-        break;
-      default: assert(0);
-      }
-      return out;
-    }
-  };
-
   //Complex float Reduce
-  template<>
-  inline Grid::ComplexF Reduce<Grid::ComplexF, u128f>::operator()(u128f in){ //2 complex
-    return Grid::ComplexF(in.f[0] + in.f[2], in.f[1] + in.f[3]);
+  template <>
+  inline Grid::ComplexF Reduce<Grid::ComplexF, vecf>::operator()(vecf in){
+    float a = 0.f, b = 0.f;
+    
+    acc(in.v, a, 0, 2, W<float>::r);
+    acc(in.v, b, 1, 2, W<float>::r);
+    
+    return Grid::ComplexF(a, b);
   }
+  
   //Real float Reduce
   template<>
-  inline Grid::RealF Reduce<Grid::RealF, u128f>::operator()(u128f in){ //4 floats
-    return in.f[0] + in.f[1] + in.f[2] + in.f[3];
+  inline Grid::RealF Reduce<Grid::RealF, vecf>::operator()(vecf in){
+    float a = 0.;
+    
+    acc(in.v, a, 0, 1, W<float>::r);
+    
+    return a;
   }
   
-  
   //Complex double Reduce
   template<>
-  inline Grid::ComplexD Reduce<Grid::ComplexD, u128d>::operator()(u128d in){ //1 complex
-    return Grid::ComplexD(in.f[0],in.f[1]);
+  inline Grid::ComplexD Reduce<Grid::ComplexD, vecd>::operator()(vecd in){
+    double a = 0., b = 0.;
+    
+    acc(in.v, a, 0, 2, W<double>::r);
+    acc(in.v, b, 1, 2, W<double>::r);
+    
+    return Grid::ComplexD(a, b);
   }
   
   //Real double Reduce
   template<>
-  inline Grid::RealD Reduce<Grid::RealD, u128d>::operator()(u128d in){ //2 doubles
-    return in.f[0] + in.f[1];
+  inline Grid::RealD Reduce<Grid::RealD, vecd>::operator()(vecd in){
+    double a = 0.f;
+    
+    acc(in.v, a, 0, 1, W<double>::r);
+    
+    return a;
   }
 
   //Integer Reduce
   template<>
   inline Integer Reduce<Integer, int>::operator()(int in){
-    // FIXME unimplemented
-   printf("Reduce : Missing integer implementation -> FIX\n");
-    assert(0);
+    return in;
   }
 }
 
 //////////////////////////////////////////////////////////////////////////////////////
 // Here assign types 
 
-  typedef Optimization::u128f SIMD_Ftype;  // Single precision type
-  typedef Optimization::u128d SIMD_Dtype; // Double precision type
+  typedef Optimization::vecf SIMD_Ftype; // Single precision type
+  typedef Optimization::vecd SIMD_Dtype; // Double precision type
   typedef int SIMD_Itype; // Integer type
 
   // prefetch utilities
   inline void v_prefetch0(int size, const char *ptr){};
   inline void prefetch_HINT_T0(const char *ptr){};
 
-
-
-  // Gpermute function
-  template < typename VectorSIMD > 
-    inline void Gpermute(VectorSIMD &y,const VectorSIMD &b, int perm ) {
-    Optimization::permute(y.v,b.v,perm);
-  }
-
-
   // Function name aliases
   typedef Optimization::Vsplat   VsplatSIMD;
   typedef Optimization::Vstore   VstoreSIMD;
@@ -481,16 +450,13 @@ namespace Optimization {
   typedef Optimization::Vstream  VstreamSIMD;
   template <typename S, typename T> using ReduceSIMD = Optimization::Reduce<S,T>;
 
- 
-
-
   // Arithmetic operations
   typedef Optimization::Sum         SumSIMD;
   typedef Optimization::Sub         SubSIMD;
+  typedef Optimization::Div         DivSIMD;
   typedef Optimization::Mult        MultSIMD;
   typedef Optimization::MultComplex MultComplexSIMD;
   typedef Optimization::Conj        ConjSIMD;
   typedef Optimization::TimesMinusI TimesMinusISIMD;
   typedef Optimization::TimesI      TimesISIMD;
-
 }

lib/simd/Grid_imci.h

@@ -26,14 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-//----------------------------------------------------------------------
-/*! @file Grid_knc.h
-  @brief Optimization libraries for AVX512 instructions set for KNC
-
-  Using intrinsics
-*/
-// Time-stamp: <2015-06-09 14:27:28 neo>
-//----------------------------------------------------------------------
 
 #include <immintrin.h>
 #include <zmmintrin.h>

lib/simd/Grid_vector_types.h

@@ -38,7 +38,7 @@ directory
 #ifndef GRID_VECTOR_TYPES
 #define GRID_VECTOR_TYPES
 
-#ifdef GENERIC_VEC
+#ifdef GEN
 #include "Grid_generic.h"
 #endif
 #ifdef SSE4

m4/ac_prog_doxygen.m4

@@ -1,54 +0,0 @@
-dnl Check for doxygen to create API docs
-dnl
-AC_DEFUN([AC_PROG_DOXYGEN],
-[
-AC_ARG_ENABLE(doxygen,
-AS_HELP_STRING([--enable-doxygen],[enable documentation generation with doxygen (auto)]))
-
-AC_ARG_ENABLE(dot,
-AS_HELP_STRING([--enable-dot],[use 'dot' to generate graphs in doxygen (auto)]))              
-AC_ARG_ENABLE(html-docs,
-AS_HELP_STRING([--enable-html-docs],[enable HTML generation with doxygen (yes)]),
-[],[ enable_html_docs=yes])              
-AC_ARG_ENABLE(latex-docs,
-AS_HELP_STRING([--enable-latex-docs],
-[enable LaTeX documentation generation with doxygen (no)]),[],[enable_latex_docs=no])              
-
-if test "x$enable_doxygen" = xno; then
- enable_doc=no	
-else 
- AC_CHECK_PROG(DOXYGEN, doxygen, doxygen)
- if test x$DOXYGEN = x; then
-  if test "x$enable_doxygen" = xyes; then
-   AC_MSG_ERROR([could not find doxygen])
-  fi
-  enable_doc=no
- else
-  doxy_ver=`doxygen --version`
-  doxy_major=`expr "$doxy_ver" : '\(@<:@0-9@:>@\)\..*'`
-  doxy_minor=`expr "$doxy_ver" : '@<:@0-9@:>@\.\(@<:@0-9@:>@\).*'`
-  if test $doxy_major -eq "1" -a $doxy_minor -ge "3" ; then
-   enable_doc=yes
-   AC_CHECK_PROG(DOT, dot, dot)
-  else
-   AC_MSG_WARN([doxygen version $doxy_ver too old, doxygen will not be used.])
-   enable_doc=no
-  fi
- fi
-fi
-
-AM_CONDITIONAL(DOXYGEN_DOC, test x$enable_doc = xyes)
-
-if test x$DOT = x; then
- if test "x$enable_dot" = xyes; then
-  AC_MSG_ERROR([could not find dot])
- fi
- enable_dot=no
-else
- enable_dot=yes
-fi
-
-AC_SUBST(enable_dot)
-AC_SUBST(enable_html_docs)
-AC_SUBST(enable_latex_docs)
-])