Merge branch 'develop' into feature/sycl

Author: jcosborn

CMakeLists.txt

@@ -192,14 +192,6 @@ option(QUDA_BUILD_ALL_TESTS "build tests by default" ON)
 option(QUDA_INSTALL_ALL_TESTS "install tests by default" ON)
 option(QUDA_BUILD_SHAREDLIB "build quda as a shared lib" ON)
 
-set(QUDA_ORDER_FP "8" CACHE STRING "which data order to use for Wilson, gauge (recon-8/9) and clover fixed-point fields (4, 8)")
-set(QUDA_ORDER_SP_MG "2" CACHE STRING "which data order to use for fp32 MG fields (2, 4)")
-set(QUDA_ORDER_FP_MG "2" CACHE STRING "which data order to use for fixed-point MG fields (2, 4, 8)")
-
-set_property(CACHE QUDA_ORDER_FP PROPERTY STRINGS 4 8)
-set_property(CACHE QUDA_ORDER_SP_MG PROPERTY STRINGS 2 4)
-set_property(CACHE QUDA_ORDER_FP_MG PROPERTY STRINGS 2 4 8)
-
 option(QUDA_BUILD_NATIVE_LAPACK "build the native blas/lapack library according to QUDA_TARGET" ON)
 option(QUDA_BUILD_NATIVE_FFT "build the native FFT library according to QUDA_TARGET" ON)
 
@@ -235,6 +227,8 @@ set(QUDA_RECONSTRUCT
   "7"
   CACHE STRING "which reconstructs to instantiate in QUDA (3-bit number - 18, 13/12, 9/8)")
 
+option(QUDA_FLUSH_DENORMALS "Whether to fliush denormals to zero" OFF)
+
 option(QUDA_CLOVER_DYNAMIC "Dynamically invert the clover term" ON)
 option(QUDA_CLOVER_RECONSTRUCT "set to ON to enable compressed clover storage (requires QUDA_CLOVER_DYNAMIC)" ON)
 option(QUDA_CLOVER_CHOLESKY_PROMOTE "Whether to promote the internal precision when inverting the clover term" ON)
@@ -248,8 +242,6 @@ option(QUDA_CTEST_DISABLE_BENCHMARKS "Disable benchmark test" ON)
 option(QUDA_FAST_COMPILE_REDUCE "enable fast compilation in blas and reduction kernels (single warp per reduction)" OFF)
 option(QUDA_FAST_COMPILE_DSLASH "enable fast compilation in dslash kernels (~20% perf impact)" OFF)
 
-option(QUDA_ALTERNATIVE_I_TO_F "enable using alternative integer-to-float conversion" OFF)
-
 option(QUDA_OPENMP "enable OpenMP" OFF)
 set(QUDA_CXX_STANDARD
     17
@@ -277,19 +269,15 @@ option(QUDA_GENERATE_DOXYGEN "generate doxygen documentation")
 mark_as_advanced(QUDA_BUILD_ALL_TESTS)
 mark_as_advanced(QUDA_INSTALL_ALL_TESTS)
 
-mark_as_advanced(QUDA_ORDER_FP)
-mark_as_advanced(QUDA_ORDER_SP_MG)
-mark_as_advanced(QUDA_ORDER_FP_MG)
 mark_as_advanced(QUDA_FAST_COMPILE_REDUCE)
 mark_as_advanced(QUDA_FAST_COMPILE_DSLASH)
 
-mark_as_advanced(QUDA_ALTERNATIVE_I_TO_F)
-
 mark_as_advanced(QUDA_MAX_MULTI_BLAS_N)
 mark_as_advanced(QUDA_MAX_MULTI_RHS)
 mark_as_advanced(QUDA_MAX_MULTI_RHS_TILE)
 mark_as_advanced(QUDA_MAX_KERNEL_ARG_SIZE)
 mark_as_advanced(QUDA_PRECISION)
+mark_as_advanced(QUDA_FLUSH_DENORMALS)
 mark_as_advanced(QUDA_RECONSTRUCT)
 mark_as_advanced(QUDA_CLOVER_CHOLESKY_PROMOTE)
 mark_as_advanced(QUDA_MULTIGRID_DSLASH_PROMOTE)

include/blas_helper.cuh

@@ -50,8 +50,8 @@ namespace quda
 
 #ifdef QUAD_SUM
   __host__ __device__ inline double set(doubledouble &a) { return a.head(); }
-  __host__ __device__ inline double2 set(doubledouble2 &a) { return make_double2(a.x.head(), a.y.head()); }
-  __host__ __device__ inline double3 set(doubledouble3 &a) { return make_double3(a.x.head(), a.y.head(), a.z.head()); }
+  __host__ __device__ inline double2 set(doubledouble2 &a) { return {a.x.head(), a.y.head()}; }
+  __host__ __device__ inline double3 set(doubledouble3 &a) { return {a.x.head(), a.y.head(), a.z.head()}; }
   __host__ __device__ inline void sum(double &a, doubledouble &b) { a += b.head(); }
   __host__ __device__ inline void sum(double2 &a, doubledouble2 &b)
   {
@@ -99,20 +99,14 @@ namespace quda
        Helper struct that contains the meta data required for
        read and writing to a spinor field in the BLAS kernels.
        @tparam store_t Type used to store field in memory
-       @tparam N Length of vector
     */
-    template <typename store_t, int N, bool is_fixed> struct data_t {
-      store_t *spinor;
-      int stride;
-      unsigned int cb_offset;
-      data_t() :
-        spinor(nullptr),
-        stride(0),
-        cb_offset(0)
-      {}
-
+    template <typename store_t, bool is_fixed> struct data_t {
+      store_t *spinor = nullptr;
+      int stride = 0;
+      unsigned int cb_offset = 0;
+      data_t() = default;
       data_t(const ColorSpinorField &x) :
-        spinor(x.data<store_t *>()), stride(x.VolumeCB()), cb_offset(x.Bytes() / (2 * sizeof(store_t) * N))
+        spinor(x.data<store_t *>()), stride(x.VolumeCB()), cb_offset(x.Bytes() / (2 * sizeof(store_t)))
       {}
     };
 
@@ -122,28 +116,20 @@ namespace quda
        specialized variant for fixed-point fields where need to store
        the meta data for the norm field.
        @tparam store_t Type used to store field in memory
-       @tparam N Length of vector
     */
-    template <typename store_t, int N> struct data_t<store_t, N, true> {
+    template <typename store_t> struct data_t<store_t, true> {
       using norm_t = float;
-      store_t *spinor;
-      norm_t *norm;
-      int stride;
-      unsigned int cb_offset;
-      unsigned int cb_norm_offset;
-      data_t() :
-        spinor(nullptr),
-        norm(nullptr),
-        stride(0),
-        cb_offset(0),
-        cb_norm_offset(0)
-      {}
-
+      store_t *spinor = nullptr;
+      norm_t *norm = nullptr;
+      int stride = 0;
+      unsigned int cb_offset = 0;
+      unsigned int cb_norm_offset = 0;
+      data_t() = default;
       data_t(const ColorSpinorField &x) :
         spinor(x.data<store_t *>()),
         norm(static_cast<norm_t *>(x.Norm())),
         stride(x.VolumeCB()),
-        cb_offset(x.Bytes() / (2 * sizeof(store_t) * N)),
+        cb_offset(x.Bytes() / (2 * sizeof(store_t))),
         cb_norm_offset(x.Bytes() / (2 * sizeof(norm_t)))
       {}
     };
@@ -156,10 +142,9 @@ namespace quda
     template <typename store_t, int N> struct Spinor {
       using Vector = typename VectorType<store_t, N>::type;
       using norm_t = float;
-      data_t<store_t, N, isFixed<store_t>::value> data;
-
-      Spinor() {}
+      data_t<store_t, isFixed<store_t>::value> data;
 
+      Spinor() = default;
       Spinor(const ColorSpinorField &x) : data(x) {}
 
       /**
@@ -234,37 +219,37 @@ namespace quda
           array<real, len> v_;
 
           constexpr int M = len / N;
+          constexpr int Nrem = len - M * N;
 #pragma unroll
           for (int i = 0; i < M; i++) {
             // first load from memory
-            Vector vecTmp = vector_load<Vector>(data.spinor, parity * data.cb_offset + x + data.stride * i);
+            auto vecTmp = vector_load<store_t, N>(data.spinor + parity * data.cb_offset, data.stride * i + x);
             // now copy into output and scale
-#pragma unroll
-            // for (int j = 0; j < N; j++) copy_and_scale(v_[i * N + j], reinterpret_cast<store_t *>(&vecTmp)[j], nrm);
-            for (int j = 0; j < N; j++) copy_and_scale(v_[i * N + j], elem(vecTmp, j), nrm);
+            copy_and_scale(&v_[i * N], vecTmp, nrm);
+          }
+          if constexpr (Nrem > 0) {
+            // first load from memory
+            auto vecTmp = vector_load<store_t, Nrem>(data.spinor + parity * data.cb_offset + data.stride * M * N, x);
+            // now copy into output and scale
+            copy_and_scale(&v_[M * N], vecTmp, nrm);
           }
 
+#pragma unroll
           for (int i = 0; i < n; i++) { v[i] = complex<real>(v_[2 * i + 0], v_[2 * i + 1]); }
         } else {
           // specialized path for half precision staggered
-          using Vector = int4;
           auto cb_offset = data.cb_norm_offset / 4;
           norm_t nrm;
           array<real, len> v_;
 
           // first load from memory
-          Vector vecTmp = vector_load<Vector>(data.spinor, parity * cb_offset + x);
+          auto vecTmp = vector_load<store_t, 8>(data.spinor, parity * cb_offset + x);
 
           // extract norm
-          memcpy(&nrm, &vecTmp.w, sizeof(norm_t));
-          array<short, 6> vecTmpShort;
-          memcpy(&vecTmpShort, &vecTmp, sizeof(vecTmpShort));
+          memcpy(&nrm, &vecTmp[6], sizeof(norm_t));
 
           // now copy into output and scale
-#pragma unroll
-          // for (int i = 0; i < len; i++) copy_and_scale(v_[i], reinterpret_cast<store_t *>(&vecTmp)[i], nrm);
-          // for (int i = 0; i < len; i++) copy_and_scale(v_[i], elem(vecTmp, i), nrm);
-          for (int i = 0; i < len; i++) copy_and_scale(v_[i], vecTmpShort[i], nrm);
+          copy_and_scale(&v_[0], vecTmp, nrm);
 
 #pragma unroll
           for (int i = 0; i < n; i++) { v[i] = complex<real>(v_[2 * i + 0], v_[2 * i + 1]); }
@@ -284,56 +269,47 @@ namespace quda
       {
         constexpr int len = 2 * n; // real-valued length
 
-        if constexpr (!(n == 3 && isHalf<store_t>::value)) {
-          array<real, len> v_;
-
-          if constexpr (isFixed<store_t>::value) {
-            real scale_inv = store_norm<isFixed<store_t>::value, real, n>(v, data.norm[x + parity * data.cb_norm_offset]);
-#pragma unroll
-            for (int i = 0; i < n; i++) {
-              v_[2 * i + 0] = scale_inv * v[i].real();
-              v_[2 * i + 1] = scale_inv * v[i].imag();
-            }
-          } else {
+        array<real, len> v_;
 #pragma unroll
-            for (int i = 0; i < n; i++) {
-              v_[2 * i + 0] = v[i].real();
-              v_[2 * i + 1] = v[i].imag();
-            }
-          }
+        for (int i = 0; i < n; i++) {
+          v_[2 * i + 0] = v[i].real();
+          v_[2 * i + 1] = v[i].imag();
+        }
+
+        if constexpr (!(n == 3 && isHalf<store_t>::value)) {
+          real scale_inv = 0.0;
+          if constexpr (isFixed<store_t>::value)
+            scale_inv = store_norm<isFixed<store_t>::value, real, n>(v, data.norm[x + parity * data.cb_norm_offset]);
 
           constexpr int M = len / N;
+          constexpr int Nrem = len - M * N;
 #pragma unroll
           for (int i = 0; i < M; i++) {
-            Vector vecTmp;
+            array<store_t, N> vecTmp;
             // first do scalar copy converting into storage type
-#pragma unroll
-            // for (int j = 0; j < N; j++) copy_scaled(reinterpret_cast<store_t *>(&vecTmp)[j], v_[i * N + j]);
-            for (int j = 0; j < N; j++) copy_scaled(elem(vecTmp, j), v_[i * N + j]);
+            copy_and_scale<store_t, real, N>(vecTmp, &v_[i * N], scale_inv);
+            // second do vectorized copy into memory
+            vector_store(data.spinor + parity * data.cb_offset, data.stride * i + x, vecTmp);
+          }
+
+          if constexpr (Nrem > 0) {
+            array<store_t, Nrem> vecTmp;
+            // first do copy converting into storage type
+            copy_and_scale<store_t, real, Nrem>(vecTmp, &v_[M * N], scale_inv);
             // second do vectorized copy into memory
-            vector_store(data.spinor, parity * data.cb_offset + x + data.stride * i, vecTmp);
+            vector_store(data.spinor + parity * data.cb_offset + data.stride * M * N, x, vecTmp);
           }
         } else {
           // specialized path for half precision staggered
-          using Vector = int4;
           auto cb_offset = data.cb_norm_offset / 4;
           norm_t norm;
           norm_t scale_inv = store_norm<isFixed<store_t>::value, real, n>(v, norm);
-          array<real, len> v_;
-#pragma unroll
-          for (int i = 0; i < n; i++) {
-            v_[2 * i + 0] = scale_inv * v[i].real();
-            v_[2 * i + 1] = scale_inv * v[i].imag();
-          }
 
-          Vector vecTmp;
-          memcpy(&vecTmp.w, &norm, sizeof(norm_t)); // pack the norm
-          array<short, 6> vecTmpShort;
-#pragma unroll
-          // for (int i = 0; i < len; i++) copy_scaled(reinterpret_cast<store_t *>(&vecTmp)[i], v_[i]);
-          // for (int i = 0; i < len; i++) copy_scaled(elem(vecTmp, i), v_[i]);
-          for (int i = 0; i < len; i++) copy_scaled(vecTmpShort[i], v_[i]);
-          memcpy(&vecTmp, &vecTmpShort, sizeof(vecTmpShort));
+          array<store_t, 8> vecTmp;
+          memcpy(&vecTmp[6], &norm, sizeof(norm_t)); // pack the norm
+          array<store_t, 6> vecTmp2;
+          copy_and_scale<store_t, real, 6>(vecTmp2, &v_[0], scale_inv);
+          std::memcpy(&vecTmp, &vecTmp2, sizeof(vecTmp2));
           // second do vectorized copy into memory
           vector_store(data.spinor, parity * cb_offset + x, vecTmp);
         }
@@ -349,42 +325,59 @@ namespace quda
        @tparam site_unroll Whether we enforce all site components must
        be unrolled onto the same thread (required for fixed-point precision)
     */
-    template <typename store_t, bool GPU, int nSpin, bool site_unroll> constexpr int n_vector() { return 0; }
+    template <typename store_t, bool GPU> constexpr int n_vector(int, int) { return 0; }
 
     // native ordering
-    template <> constexpr int n_vector<double, true, 4, false>() { return 2; }
-    template <> constexpr int n_vector<double, true, 1, false>() { return 2; }
-
-    template <> constexpr int n_vector<double, true, 4, true>() { return 2; }
-    template <> constexpr int n_vector<double, true, 1, true>() { return 2; }
-
-    template <> constexpr int n_vector<float, true, 4, false>() { return 4; }
-    template <> constexpr int n_vector<float, true, 1, false>() { return 4; }
+    template <> constexpr int n_vector<double, true>(int nSpin, int site_unroll)
+    {
+      if (site_unroll)
+        return nSpin == 4 ? colorspinor::get_vector_order<double>(24) : colorspinor::get_vector_order<double>(6);
+      else
+        return colorspinor::get_vector_order<double>(4);
+    }
 
-    template <> constexpr int n_vector<float, true, 4, true>() { return 4; }
-    template <> constexpr int n_vector<float, true, 1, true>() { return 2; }
+    template <> constexpr int n_vector<float, true>(int nSpin, int site_unroll)
+    {
+      if (site_unroll)
+        return nSpin == 4 ? colorspinor::get_vector_order<float>(24) : colorspinor::get_vector_order<float>(6);
+      else
+        return colorspinor::get_vector_order<float>(8);
+    }
 
-    template <> constexpr int n_vector<short, true, 4, true>() { return QUDA_ORDER_FP; }
-    template <> constexpr int n_vector<short, true, 1, true>() { return 2; }
+    template <> constexpr int n_vector<short, true>(int nSpin, int site_unroll)
+    {
+      if (site_unroll)
+        return nSpin == 4 ? colorspinor::get_vector_order<short>(24) : colorspinor::get_vector_order<short>(6);
+      else
+        return colorspinor::get_vector_order<short>(16);
+    }
 
-    template <> constexpr int n_vector<int8_t, true, 4, true>() { return QUDA_ORDER_FP; }
-    template <> constexpr int n_vector<int8_t, true, 1, true>() { return 2; }
+    template <> constexpr int n_vector<int8_t, true>(int nSpin, int site_unroll)
+    {
+      if (site_unroll)
+        return nSpin == 4 ? colorspinor::get_vector_order<int8_t>(24) : colorspinor::get_vector_order<int8_t>(6);
+      else
+        return colorspinor::get_vector_order<int8_t>(16);
+    }
 
     // Just use float-2/float-4 ordering on CPU when not site unrolling
-    template <> constexpr int n_vector<double, false, 4, false>() { return 2; }
-    template <> constexpr int n_vector<double, false, 1, false>() { return 2; }
-    template <> constexpr int n_vector<float, false, 4, false>() { return 4; }
-    template <> constexpr int n_vector<float, false, 1, false>() { return 4; }
-
-    // AoS ordering is used on CPU uses when we are site unrolling
-    template <> constexpr int n_vector<double, false, 4, true>() { return 24; }
-    template <> constexpr int n_vector<double, false, 1, true>() { return 6; }
-    template <> constexpr int n_vector<float, false, 4, true>() { return 24; }
-    template <> constexpr int n_vector<float, false, 1, true>() { return 6; }
-    template <> constexpr int n_vector<short, false, 4, true>() { return 24; }
-    template <> constexpr int n_vector<short, false, 1, true>() { return 6; }
-    template <> constexpr int n_vector<int8_t, false, 4, true>() { return 24; }
-    template <> constexpr int n_vector<int8_t, false, 1, true>() { return 6; }
+    template <> constexpr int n_vector<double, false>(int nSpin, int site_unroll)
+    {
+      if (site_unroll) {
+        return nSpin * 6;
+      } else {
+        return 2;
+      }
+    }
+
+    template <> constexpr int n_vector<float, false>(int nSpin, int site_unroll)
+    {
+      if (site_unroll) {
+        return nSpin * 6;
+      } else {
+        return 4;
+      }
+    }
 
     template <template <typename...> class Functor,
               template <template <typename...> class, typename store_t, typename y_store_t, int, typename> class Blas,