]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Specialize EvaluatorTensorProduct<evaluate_evenodd> for non-templated execution 14237/head
authorPeter Munch <peterrmuench@gmail.com>
Sat, 3 Sep 2022 12:52:50 +0000 (14:52 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Sat, 3 Sep 2022 13:43:19 +0000 (15:43 +0200)
include/deal.II/matrix_free/tensor_product_kernels.h

index a83591563dbe121c3ea28a5a9bcf30ee9b66b1e4..10971a272fd59c36d5dcee163eeba90b843cb598 100644 (file)
@@ -1567,6 +1567,218 @@ namespace internal
 
 
 
+  template <int dim,
+            int n_rows_static,
+            int n_columns_static,
+            typename Number,
+            typename Number2,
+            int  direction,
+            bool contract_over_rows,
+            bool add,
+            int  type,
+            bool one_line>
+  inline void
+  even_odd_apply(const int                       n_rows_in,
+                 const int                       n_columns_in,
+                 const Number2 *DEAL_II_RESTRICT shapes,
+                 const Number *                  in,
+                 Number *                        out)
+  {
+    static_assert(type < 3, "Only three variants type=0,1,2 implemented");
+    static_assert(one_line == false || direction == dim - 1,
+                  "Single-line evaluation only works for direction=dim-1.");
+
+    const int n_rows = n_rows_static == -1 ? n_rows_in : n_rows_static;
+    const int n_columns =
+      n_columns_static == -1 ? n_columns_in : n_columns_static;
+
+    Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
+             in != out,
+           ExcMessage("In-place operation only supported for "
+                      "n_rows==n_columns or single-line interpolation"));
+
+    // We cannot statically assert that direction is less than dim, so must do
+    // an additional dynamic check
+    AssertIndexRange(direction, dim);
+
+    const int     nn = contract_over_rows ? n_columns : n_rows;
+    const int     mm = contract_over_rows ? n_rows : n_columns;
+    constexpr int mm_static =
+      contract_over_rows ? n_rows_static : n_columns_static;
+    const int     n_cols     = nn / 2;
+    const int     mid        = mm / 2;
+    constexpr int mid_static = mm_static / 2;
+    constexpr int max_mid    = 15; // for non-templated execution
+
+    Assert((n_rows_static != -1 && n_columns_static != -1) || mid <= max_mid,
+           ExcNotImplemented());
+
+    const int stride    = Utilities::pow(n_columns, direction);
+    const int n_blocks1 = one_line ? 1 : stride;
+    const int n_blocks2 =
+      Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
+
+    const int offset = (n_columns + 1) / 2;
+
+    // this code may look very inefficient at first sight due to the many
+    // different cases with if's at the innermost loop part, but all of the
+    // conditionals can be evaluated at compile time because they are
+    // templates, so the compiler should optimize everything away
+    for (int i2 = 0; i2 < n_blocks2; ++i2)
+      {
+        for (int i1 = 0; i1 < n_blocks1; ++i1)
+          {
+            constexpr unsigned int mid_size =
+              (n_rows_static == -1 || n_columns_static == -1) ?
+                max_mid :
+                (mid_static > 0 ? mid_static : 1);
+            Number xp[mid_size], xm[mid_size];
+            for (int i = 0; i < mid; ++i)
+              {
+                if (contract_over_rows == true && type == 1)
+                  {
+                    xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
+                    xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
+                  }
+                else
+                  {
+                    xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
+                    xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
+                  }
+              }
+            Number xmid = in[stride * mid];
+            for (int col = 0; col < n_cols; ++col)
+              {
+                Number r0, r1;
+                if (mid > 0)
+                  {
+                    if (contract_over_rows == true)
+                      {
+                        r0 = shapes[col] * xp[0];
+                        r1 = shapes[(n_rows - 1) * offset + col] * xm[0];
+                      }
+                    else
+                      {
+                        r0 = shapes[col * offset] * xp[0];
+                        r1 = shapes[(n_rows - 1 - col) * offset] * xm[0];
+                      }
+                    for (int ind = 1; ind < mid; ++ind)
+                      {
+                        if (contract_over_rows == true)
+                          {
+                            r0 += shapes[ind * offset + col] * xp[ind];
+                            r1 += shapes[(n_rows - 1 - ind) * offset + col] *
+                                  xm[ind];
+                          }
+                        else
+                          {
+                            r0 += shapes[col * offset + ind] * xp[ind];
+                            r1 += shapes[(n_rows - 1 - col) * offset + ind] *
+                                  xm[ind];
+                          }
+                      }
+                  }
+                else
+                  r0 = r1 = Number();
+                if (mm % 2 == 1 && contract_over_rows == true)
+                  {
+                    if (type == 1)
+                      r1 += shapes[mid * offset + col] * xmid;
+                    else
+                      r0 += shapes[mid * offset + col] * xmid;
+                  }
+                else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0 || mm == 3))
+                  r0 += shapes[col * offset + mid] * xmid;
+
+                if (add)
+                  {
+                    out[stride * col] += r0 + r1;
+                    if (type == 1 && contract_over_rows == false)
+                      out[stride * (nn - 1 - col)] += r1 - r0;
+                    else
+                      out[stride * (nn - 1 - col)] += r0 - r1;
+                  }
+                else
+                  {
+                    out[stride * col] = r0 + r1;
+                    if (type == 1 && contract_over_rows == false)
+                      out[stride * (nn - 1 - col)] = r1 - r0;
+                    else
+                      out[stride * (nn - 1 - col)] = r0 - r1;
+                  }
+              }
+            if (type == 0 && contract_over_rows == true && nn % 2 == 1 &&
+                mm % 2 == 1 && mm > 3)
+              {
+                if (add)
+                  out[stride * n_cols] += shapes[mid * offset + n_cols] * xmid;
+                else
+                  out[stride * n_cols] = shapes[mid * offset + n_cols] * xmid;
+              }
+            else if (contract_over_rows == true && nn % 2 == 1)
+              {
+                Number r0;
+                if (mid > 0)
+                  {
+                    r0 = shapes[n_cols] * xp[0];
+                    for (int ind = 1; ind < mid; ++ind)
+                      r0 += shapes[ind * offset + n_cols] * xp[ind];
+                  }
+                else
+                  r0 = Number();
+                if (type != 1 && mm % 2 == 1)
+                  r0 += shapes[mid * offset + n_cols] * xmid;
+
+                if (add)
+                  out[stride * n_cols] += r0;
+                else
+                  out[stride * n_cols] = r0;
+              }
+            else if (contract_over_rows == false && nn % 2 == 1)
+              {
+                Number r0;
+                if (mid > 0)
+                  {
+                    if (type == 1)
+                      {
+                        r0 = shapes[n_cols * offset] * xm[0];
+                        for (int ind = 1; ind < mid; ++ind)
+                          r0 += shapes[n_cols * offset + ind] * xm[ind];
+                      }
+                    else
+                      {
+                        r0 = shapes[n_cols * offset] * xp[0];
+                        for (int ind = 1; ind < mid; ++ind)
+                          r0 += shapes[n_cols * offset + ind] * xp[ind];
+                      }
+                  }
+                else
+                  r0 = Number();
+
+                if ((type == 0 || type == 2) && mm % 2 == 1)
+                  r0 += shapes[n_cols * offset + mid] * xmid;
+
+                if (add)
+                  out[stride * n_cols] += r0;
+                else
+                  out[stride * n_cols] = r0;
+              }
+            if (one_line == false)
+              {
+                in += 1;
+                out += 1;
+              }
+          }
+        if (one_line == false)
+          {
+            in += stride * (mm - 1);
+            out += stride * (nn - 1);
+          }
+      }
+  }
+
+
+
   /**
    * Internal evaluator for 1d-3d shape function using the tensor product form
    * of the basis functions.
@@ -1751,7 +1963,19 @@ namespace internal
     static void
     apply(const Number2 *DEAL_II_RESTRICT shape_data,
           const Number *                  in,
-          Number *                        out);
+          Number *                        out)
+    {
+      even_odd_apply<dim,
+                     n_rows,
+                     n_columns,
+                     Number,
+                     Number2,
+                     direction,
+                     contract_over_rows,
+                     add,
+                     type,
+                     one_line>(n_rows, n_columns, shape_data, in, out);
+    }
 
   private:
     const Number2 *shape_values;
@@ -1760,203 +1984,133 @@ namespace internal
   };
 
 
-
-  template <int dim,
-            int n_rows,
-            int n_columns,
-            typename Number,
-            typename Number2>
-  template <int  direction,
-            bool contract_over_rows,
-            bool add,
-            int  type,
-            bool one_line>
-  inline void
-  EvaluatorTensorProduct<evaluate_evenodd,
-                         dim,
-                         n_rows,
-                         n_columns,
-                         Number,
-                         Number2>::apply(const Number2 *DEAL_II_RESTRICT shapes,
-                                         const Number *                  in,
-                                         Number *                        out)
+  /**
+   * Internal evaluator for shape function using the tensor product form
+   * of the basis functions. The same as the other templated class but
+   * without making use of template arguments and variable loop bounds
+   * instead.
+   */
+  template <int dim, typename Number, typename Number2>
+  struct EvaluatorTensorProduct<evaluate_evenodd, dim, 0, 0, Number, Number2>
   {
-    static_assert(type < 3, "Only three variants type=0,1,2 implemented");
-    static_assert(one_line == false || direction == dim - 1,
-                  "Single-line evaluation only works for direction=dim-1.");
-    Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
-             in != out,
-           ExcMessage("In-place operation only supported for "
-                      "n_rows==n_columns or single-line interpolation"));
+    EvaluatorTensorProduct()
+      : shape_values(nullptr)
+      , shape_gradients(nullptr)
+      , shape_hessians(nullptr)
+      , n_rows(numbers::invalid_unsigned_int)
+      , n_columns(numbers::invalid_unsigned_int)
+    {}
 
-    // We cannot statically assert that direction is less than dim, so must do
-    // an additional dynamic check
-    AssertIndexRange(direction, dim);
+    EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
+                           const unsigned int            n_rows    = 0,
+                           const unsigned int            n_columns = 0)
+      : shape_values(shape_values.begin())
+      , shape_gradients(nullptr)
+      , shape_hessians(nullptr)
+    {
+      AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
+    }
 
-    constexpr int nn     = contract_over_rows ? n_columns : n_rows;
-    constexpr int mm     = contract_over_rows ? n_rows : n_columns;
-    constexpr int n_cols = nn / 2;
-    constexpr int mid    = mm / 2;
+    EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
+                           const AlignedVector<Number2> &shape_gradients,
+                           const AlignedVector<Number2> &shape_hessians,
+                           const unsigned int            n_rows    = 0,
+                           const unsigned int            n_columns = 0)
+      : shape_values(shape_values.begin())
+      , shape_gradients(shape_gradients.begin())
+      , shape_hessians(shape_hessians.begin())
+      , n_rows(n_rows)
+      , n_columns(n_columns)
+    {
+      if (!shape_values.empty())
+        AssertDimension(shape_values.size(), n_rows * ((n_columns + 1) / 2));
+      if (!shape_gradients.empty())
+        AssertDimension(shape_gradients.size(), n_rows * ((n_columns + 1) / 2));
+      if (!shape_hessians.empty())
+        AssertDimension(shape_hessians.size(), n_rows * ((n_columns + 1) / 2));
+    }
 
-    constexpr int stride    = Utilities::pow(n_columns, direction);
-    constexpr int n_blocks1 = one_line ? 1 : stride;
-    constexpr int n_blocks2 =
-      Utilities::pow(n_rows, (direction >= dim) ? 0 : (dim - direction - 1));
+    template <int direction, bool contract_over_rows, bool add>
+    void
+    values(const Number in[], Number out[]) const
+    {
+      Assert(shape_values != nullptr, ExcNotInitialized());
+      apply<direction, contract_over_rows, add, 0>(shape_values, in, out);
+    }
 
-    constexpr int offset = (n_columns + 1) / 2;
+    template <int direction, bool contract_over_rows, bool add>
+    void
+    gradients(const Number in[], Number out[]) const
+    {
+      Assert(shape_gradients != nullptr, ExcNotInitialized());
+      apply<direction, contract_over_rows, add, 1>(shape_gradients, in, out);
+    }
 
-    // this code may look very inefficient at first sight due to the many
-    // different cases with if's at the innermost loop part, but all of the
-    // conditionals can be evaluated at compile time because they are
-    // templates, so the compiler should optimize everything away
-    for (int i2 = 0; i2 < n_blocks2; ++i2)
-      {
-        for (int i1 = 0; i1 < n_blocks1; ++i1)
-          {
-            Number xp[mid > 0 ? mid : 1], xm[mid > 0 ? mid : 1];
-            for (int i = 0; i < mid; ++i)
-              {
-                if (contract_over_rows == true && type == 1)
-                  {
-                    xp[i] = in[stride * i] - in[stride * (mm - 1 - i)];
-                    xm[i] = in[stride * i] + in[stride * (mm - 1 - i)];
-                  }
-                else
-                  {
-                    xp[i] = in[stride * i] + in[stride * (mm - 1 - i)];
-                    xm[i] = in[stride * i] - in[stride * (mm - 1 - i)];
-                  }
-              }
-            Number xmid = in[stride * mid];
-            for (int col = 0; col < n_cols; ++col)
-              {
-                Number r0, r1;
-                if (mid > 0)
-                  {
-                    if (contract_over_rows == true)
-                      {
-                        r0 = shapes[col] * xp[0];
-                        r1 = shapes[(n_rows - 1) * offset + col] * xm[0];
-                      }
-                    else
-                      {
-                        r0 = shapes[col * offset] * xp[0];
-                        r1 = shapes[(n_rows - 1 - col) * offset] * xm[0];
-                      }
-                    for (int ind = 1; ind < mid; ++ind)
-                      {
-                        if (contract_over_rows == true)
-                          {
-                            r0 += shapes[ind * offset + col] * xp[ind];
-                            r1 += shapes[(n_rows - 1 - ind) * offset + col] *
-                                  xm[ind];
-                          }
-                        else
-                          {
-                            r0 += shapes[col * offset + ind] * xp[ind];
-                            r1 += shapes[(n_rows - 1 - col) * offset + ind] *
-                                  xm[ind];
-                          }
-                      }
-                  }
-                else
-                  r0 = r1 = Number();
-                if (mm % 2 == 1 && contract_over_rows == true)
-                  {
-                    if (type == 1)
-                      r1 += shapes[mid * offset + col] * xmid;
-                    else
-                      r0 += shapes[mid * offset + col] * xmid;
-                  }
-                else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0 || mm == 3))
-                  r0 += shapes[col * offset + mid] * xmid;
+    template <int direction, bool contract_over_rows, bool add>
+    void
+    hessians(const Number in[], Number out[]) const
+    {
+      Assert(shape_hessians != nullptr, ExcNotInitialized());
+      apply<direction, contract_over_rows, add, 2>(shape_hessians, in, out);
+    }
 
-                if (add)
-                  {
-                    out[stride * col] += r0 + r1;
-                    if (type == 1 && contract_over_rows == false)
-                      out[stride * (nn - 1 - col)] += r1 - r0;
-                    else
-                      out[stride * (nn - 1 - col)] += r0 - r1;
-                  }
-                else
-                  {
-                    out[stride * col] = r0 + r1;
-                    if (type == 1 && contract_over_rows == false)
-                      out[stride * (nn - 1 - col)] = r1 - r0;
-                    else
-                      out[stride * (nn - 1 - col)] = r0 - r1;
-                  }
-              }
-            if (type == 0 && contract_over_rows == true && nn % 2 == 1 &&
-                mm % 2 == 1 && mm > 3)
-              {
-                if (add)
-                  out[stride * n_cols] += shapes[mid * offset + n_cols] * xmid;
-                else
-                  out[stride * n_cols] = shapes[mid * offset + n_cols] * xmid;
-              }
-            else if (contract_over_rows == true && nn % 2 == 1)
-              {
-                Number r0;
-                if (mid > 0)
-                  {
-                    r0 = shapes[n_cols] * xp[0];
-                    for (int ind = 1; ind < mid; ++ind)
-                      r0 += shapes[ind * offset + n_cols] * xp[ind];
-                  }
-                else
-                  r0 = Number();
-                if (type != 1 && mm % 2 == 1)
-                  r0 += shapes[mid * offset + n_cols] * xmid;
+    template <int direction, bool contract_over_rows, bool add>
+    void
+    values_one_line(const Number in[], Number out[]) const
+    {
+      Assert(shape_values != nullptr, ExcNotInitialized());
+      apply<direction, contract_over_rows, add, 0, true>(shape_values, in, out);
+    }
 
-                if (add)
-                  out[stride * n_cols] += r0;
-                else
-                  out[stride * n_cols] = r0;
-              }
-            else if (contract_over_rows == false && nn % 2 == 1)
-              {
-                Number r0;
-                if (mid > 0)
-                  {
-                    if (type == 1)
-                      {
-                        r0 = shapes[n_cols * offset] * xm[0];
-                        for (int ind = 1; ind < mid; ++ind)
-                          r0 += shapes[n_cols * offset + ind] * xm[ind];
-                      }
-                    else
-                      {
-                        r0 = shapes[n_cols * offset] * xp[0];
-                        for (int ind = 1; ind < mid; ++ind)
-                          r0 += shapes[n_cols * offset + ind] * xp[ind];
-                      }
-                  }
-                else
-                  r0 = Number();
+    template <int direction, bool contract_over_rows, bool add>
+    void
+    gradients_one_line(const Number in[], Number out[]) const
+    {
+      Assert(shape_gradients != nullptr, ExcNotInitialized());
+      apply<direction, contract_over_rows, add, 1, true>(shape_gradients,
+                                                         in,
+                                                         out);
+    }
 
-                if ((type == 0 || type == 2) && mm % 2 == 1)
-                  r0 += shapes[n_cols * offset + mid] * xmid;
+    template <int direction, bool contract_over_rows, bool add>
+    void
+    hessians_one_line(const Number in[], Number out[]) const
+    {
+      Assert(shape_hessians != nullptr, ExcNotInitialized());
+      apply<direction, contract_over_rows, add, 2, true>(shape_hessians,
+                                                         in,
+                                                         out);
+    }
 
-                if (add)
-                  out[stride * n_cols] += r0;
-                else
-                  out[stride * n_cols] = r0;
-              }
-            if (one_line == false)
-              {
-                in += 1;
-                out += 1;
-              }
-          }
-        if (one_line == false)
-          {
-            in += stride * (mm - 1);
-            out += stride * (nn - 1);
-          }
-      }
-  }
+    template <int  direction,
+              bool contract_over_rows,
+              bool add,
+              int  type,
+              bool one_line = false>
+    void
+    apply(const Number2 *DEAL_II_RESTRICT shape_data,
+          const Number *                  in,
+          Number *                        out) const
+    {
+      even_odd_apply<dim,
+                     -1,
+                     -1,
+                     Number,
+                     Number2,
+                     direction,
+                     contract_over_rows,
+                     add,
+                     type,
+                     one_line>(n_rows, n_columns, shape_data, in, out);
+    }
+
+  private:
+    const Number2 *    shape_values;
+    const Number2 *    shape_gradients;
+    const Number2 *    shape_hessians;
+    const unsigned int n_rows;
+    const unsigned int n_columns;
+  };
 
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.