Vectorized (SIMD) Numerical Schemes

Common/include/geometry/dual_grid/CEdge.hpp

@@ -35,12 +35,13 @@
  * \author F. Palacios
  */
 class CEdge {
-  static_assert(su2activematrix::Storage == StorageType::RowMajor, "Needed to return normal as pointer.");
-
+  static_assert(su2activematrix::IsRowMajor, "Needed to return normal as pointer.");
 private:
-  su2matrix<unsigned long> Nodes; /*!< \brief Vector to store the node indices of the edge. */
-  su2activematrix Normal;         /*!< \brief Normal (area) of the edge. */
-  su2activematrix Coord_CG;       /*!< \brief Center-of-gravity (mid point) of the edge. */
+  using Index = unsigned long;
+  using NodeArray = C2DContainer<Index, Index, StorageType::ColumnMajor, 64, DynamicSize, 2>;
+  NodeArray Nodes;           /*!< \brief Vector to store the node indices of the edge. */
+  su2activematrix Normal;    /*!< \brief Normal (area) of the edge. */
+  su2activematrix Coord_CG;  /*!< \brief Center-of-gravity (mid point) of the edge. */
 
 public:
   enum NodePosition : unsigned long {LEFT = 0, RIGHT = 1};
@@ -84,6 +85,14 @@ class CEdge {
    */
   inline unsigned long GetNode(unsigned long iEdge, unsigned long iNode) const { return Nodes(iEdge,iNode); }
 
+  /*!
+   * \brief SIMD version of GetNode, iNode returned for multiple sequential iEdges
+   */
+  template<class IndexSIMD_t>
+  FORCEINLINE IndexSIMD_t GetNode(unsigned long iEdge, unsigned long iNode) const {
+    return IndexSIMD_t(&Nodes(iEdge,iNode));
+  }
+
   /*!
    * \brief Set the node indices of an edge.
    * \param[in] iEdge - Edge index.
@@ -168,6 +177,11 @@ class CEdge {
    */
   inline const su2double* GetNormal(unsigned long iEdge) const { return Normal[iEdge]; }
 
+  /*!
+   * \brief Get the entire matrix of edge normals.
+   */
+  inline const su2activematrix& GetNormal() const { return Normal; }
+
   /*!
    * \brief Initialize normal vector to 0.
    */

Common/include/geometry/dual_grid/CPoint.hpp

@@ -136,6 +136,11 @@ class CPoint {
    */
   inline su2double *GetCoord(unsigned long iPoint) { return Coord[iPoint]; }
 
+  /*!
+   * \brief Get the entire matrix of coordinates of the control volumes.
+   */
+  inline const su2activematrix& GetCoord() const { return Coord; }
+
   /*!
    * \brief Set the coordinates for the control volume.
    * \param[in] iPoint - Index of the point.

Common/include/interface_interpolation/CRadialBasisFunction.hpp

@@ -34,8 +34,7 @@
  * \brief Radial basis function interpolation.
  */
 class CRadialBasisFunction final : public CInterpolator {
-  static_assert(su2passivematrix::Storage == StorageType::RowMajor,
-                "This class relies on row major storage throughout.");
+  static_assert(su2passivematrix::IsRowMajor, "This class relies on row major storage throughout.");
 private:
   unsigned long MinDonors = 0, AvgDonors = 0, MaxDonors = 0;
   passivedouble Density = 0.0, AvgCorrection = 0.0, MaxCorrection = 0.0;

Common/include/linear_algebra/CSysMatrix.hpp

@@ -30,6 +30,7 @@
 
 #include "../../include/mpi_structure.hpp"
 #include "../../include/omp_structure.hpp"
+#include "../../include/parallelization/vectorization.hpp"
 #include "CSysVector.hpp"
 #include "CPastixWrapper.hpp"
 
@@ -558,6 +559,40 @@ class CSysMatrix {
     UpdateBlocks<OtherType,-1>(iEdge, iPoint, jPoint, block_i, block_j);
   }
 
+  /*!
+   * \brief SIMD version, does the update for multiple edges and points.
+   */
+  template<class MatrixSIMD_t, class T, size_t N>
+  FORCEINLINE void UpdateBlocks(simd::Array<T,N> iEdge, simd::Array<T,N> iPoint, simd::Array<T,N> jPoint,
+                                const MatrixSIMD_t& block_i, const MatrixSIMD_t& block_j) {
+
+    /*--- Fetch the blocks for all edges. ---*/
+    ScalarType* bii[N] = {nullptr};
+    ScalarType* bjj[N] = {nullptr};
+    ScalarType* bij[N] = {nullptr};
+    ScalarType* bji[N] = {nullptr};
+
+    for (auto k = 0ul; k < N; ++k) {
+      bii[k] = &matrix[dia_ptr[iPoint[k]]*block_i.size()];
+      bjj[k] = &matrix[dia_ptr[jPoint[k]]*block_i.size()];
+      bij[k] = &matrix[edge_ptr(iEdge[k],0)*block_i.size()];
+      bji[k] = &matrix[edge_ptr(iEdge[k],1)*block_i.size()];
+    }
+
+    /*--- Unpack the SIMD elements of the input blocks. ---*/
+    for (auto iVar = 0ul; iVar < block_i.rows(); ++iVar) {
+      for (auto jVar = 0ul; jVar < block_i.cols(); ++jVar) {
+        SU2_OMP_SIMD
+        for (auto k = 0ul; k < N; ++k) {
+          *(bii[k]++) += PassiveAssign(block_i(iVar,jVar)[k]);
+          *(bij[k]++) =  PassiveAssign(block_j(iVar,jVar)[k]);
+          *(bji[k]++) = -PassiveAssign(block_i(iVar,jVar)[k]);
+          *(bjj[k]++) -= PassiveAssign(block_j(iVar,jVar)[k]);
+        }
+      }
+    }
+  }
+
   /*!
    * \brief Sets 2 blocks ij and ji (add to i* sub from j*) associated with
    *        one edge of an FVM-type sparse pattern.
@@ -602,6 +637,33 @@ class CSysMatrix {
     SetBlocks<OtherType,-1,false>(iEdge, block_i, block_j);
   }
 
+  /*!
+   * \brief SIMD version, does the update for multiple edges.
+   */
+  template<class MatrixSIMD_t, class T, size_t N>
+  FORCEINLINE void SetBlocks(simd::Array<T,N> iEdge, const MatrixSIMD_t& block_i, const MatrixSIMD_t& block_j) {
+
+    /*--- Fetch blocks for all edges. ---*/
+    ScalarType* bij[N] = {nullptr};
+    ScalarType* bji[N] = {nullptr};
+
+    for (auto k = 0ul; k < N; ++k) {
+      bij[k] = &matrix[edge_ptr(iEdge[k],0)*block_i.size()];
+      bji[k] = &matrix[edge_ptr(iEdge[k],1)*block_i.size()];
+    }
+
+    /*--- Unpack the SIMD elements of the input blocks. ---*/
+    for (auto iVar = 0ul; iVar < block_i.rows(); ++iVar) {
+      for (auto jVar = 0ul; jVar < block_i.cols(); ++jVar) {
+        SU2_OMP_SIMD
+        for (auto k = 0ul; k < N; ++k) {
+          *(bij[k]++) = PassiveAssign(block_j(iVar,jVar)[k]);
+          *(bji[k]++) = -PassiveAssign(block_i(iVar,jVar)[k]);
+        }
+      }
+    }
+  }
+
   /*!
    * \brief Sets the specified block to the (i, i) subblock of the sparse matrix.
    *        Scales the input block by factor alpha. If the Overwrite parameter is

Common/include/linear_algebra/CSysVector.hpp

@@ -31,6 +31,9 @@
 #include <cmath>
 #include <cstdlib>
 
+#include "../basic_types/datatype_structure.hpp"
+#include "../omp_structure.hpp"
+#include "../parallelization/vectorization.hpp"
 
 /*!
  * \class CSysVector
@@ -329,6 +332,41 @@ class CSysVector {
       vec_val[val_ipoint*nVar+iVar] = val_residual[iVar];
   }
 
+  /*!
+   * \brief Vectorized version of SetBlock, sets multiple iPoint's.
+   * \param[in] Overwrite - True: write over existing data; False: add to existing data.
+   * \param[in] iPoint - SIMD integer, the positions to update.
+   * \param[in] vector - Vector of SIMD scalars.
+   * \param[in] alpha - Optional scale factor (axpy type operation).
+   */
+  template<class T, size_t N, class VectorSIMD_t, bool Overwrite = true>
+  FORCEINLINE void SetBlock(simd::Array<T,N> iPoint, const VectorSIMD_t& vector, ScalarType alpha = 1) {
+    const auto nVar = vector.rows();
+    for (auto iVar = 0ul; iVar < nVar; ++iVar) {
+      SU2_OMP_SIMD
+      for (auto k = 0ul; k < N; ++k) {
+        vec_val[iPoint[k]*nVar+iVar] *= 1-Overwrite;
+        vec_val[iPoint[k]*nVar+iVar] += alpha * vector(iVar)[k];
+      }
+    }
+  }
+
+  /*!
+   * \brief Vectorized version of AddBlock, see SetBlock.
+   */
+  template<class T, size_t N, class VectorSIMD_t>
+  FORCEINLINE void AddBlock(simd::Array<T,N> iPoint, const VectorSIMD_t& vector, ScalarType alpha = 1) {
+    SetBlock<T, N, VectorSIMD_t, false>(iPoint, vector, alpha);
+  }
+
+  /*!
+   * \brief Vectorized version of SubtractBlock, see SetBlock.
+   */
+  template<class T, size_t N, class VectorSIMD_t>
+  FORCEINLINE void SubtractBlock(simd::Array<T,N> iPoint, const VectorSIMD_t& vector) {
+    SetBlock<T, N, VectorSIMD_t, false>(iPoint, vector, -1);
+  }
+
   /*!
    * \brief Set the residual to zero.
    * \param[in] val_ipoint - index of the point where set the residual.