From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Fri, 8 Sep 2017 07:02:25 +0000 (+0200)
Subject: Rewrite TensorProductPolynomials::compute for more performance
X-Git-Tag: v9.0.0-rc1~1101^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F5047%2Fhead;p=dealii.git

Rewrite TensorProductPolynomials::compute for more performance
---

diff --git a/source/base/tensor_product_polynomials.cc b/source/base/tensor_product_polynomials.cc
index 8c7a96c3d1..f7bd60aa12 100644
--- a/source/base/tensor_product_polynomials.cc
+++ b/source/base/tensor_product_polynomials.cc
@@ -18,6 +18,7 @@
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/table.h>
 
+#include <boost/container/small_vector.hpp>
 #include <array>
 
 DEAL_II_NAMESPACE_OPEN
@@ -245,6 +246,7 @@ compute (const Point<dim>            &p,
          std::vector<Tensor<3,dim> > &third_derivatives,
          std::vector<Tensor<4,dim> > &fourth_derivatives) const
 {
+  Assert(dim<=3, ExcNotImplemented());
   Assert (values.size()==n_tensor_pols    || values.size()==0,
           ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
   Assert (grads.size()==n_tensor_pols     || grads.size()==0,
@@ -262,9 +264,7 @@ compute (const Point<dim>            &p,
              update_3rd_derivatives      = (third_derivatives.size()==n_tensor_pols),
              update_4th_derivatives = (fourth_derivatives.size()==n_tensor_pols);
 
-  // check how many
-  // values/derivatives we have to
-  // compute
+  // check how many values/derivatives we have to compute
   unsigned int n_values_and_derivatives = 0;
   if (update_values)
     n_values_and_derivatives = 1;
@@ -277,159 +277,110 @@ compute (const Point<dim>            &p,
   if (update_4th_derivatives)
     n_values_and_derivatives = 5;
 
-  // Provide a shortcut if only values are requested. For this case usually the
-  // temporary memory allocation below does not pay off. Also we loop over all
-  // tensor polynomials in a peculiar way to avoid all dynamic memory allocation.
-  // We need to evaluate the polynomial value n_polynomials*dim times, and
-  // need to compute n_polynomials^dim tensor polynomials. Therefore we can split
-  // the loop over the tensor polynomials into one loop over n_polynomials,
-  // evaluate this polynomial for each dimension and multiply it with the
-  // associated n_polynomials^(dim-1) tensor polynomials before moving to the
-  // next polynomial.
-  if (n_values_and_derivatives == 1)
-    {
-      const unsigned int n_polynomials = polynomials.size();
-      const unsigned int factor = n_tensor_pols/n_polynomials;
-
-      std::fill(values.begin(),values.end(),1.0);
-
-      for (unsigned int polynomial=0; polynomial<n_polynomials; ++polynomial)
-        {
-          for (unsigned int d=0; d<dim; ++d)
-            {
-              const double value = polynomials[polynomial].value(p(d));
-
-              for (unsigned int f=0; f<factor; ++f)
-                {
-                  unsigned int index = 0;
-                  switch (d)
-                    {
-                    case 0:
-                      index = polynomial+f*n_polynomials;
-                      break;
-                    case 1:
-                      index = polynomial*n_polynomials + (f/n_polynomials)*n_polynomials*n_polynomials + f%n_polynomials;
-                      break;
-                    case 2:
-                      index = polynomial*factor + f;
-                      break;
-                    default:
-                      AssertThrow(false, ExcNotImplemented());
-                    }
-
-                  const unsigned int i = index_map_inverse[index];
-                  values[i] *= value;
-                }
-            }
-        }
-      return;
-    }
-
-
-  // Compute the values (and derivatives, if
-  // necessary) of all polynomials at this
-  // evaluation point. To avoid expensive memory allocation,
-  // we use a small amount of memory on the stack, and store the
-  // result in an array (that has enough
-  // fields for any evaluation of values and
-  // derivatives, up to the 4th derivative, for up to 20 polynomials).
-  // If someone uses a larger number of
-  // polynomials, we need to allocate more memory on the heap.
-  std::array<std::array<double,5>, dim> *v;
-  std::array<std::array<std::array<double,5>, dim>, 20> small_array;
-  std::vector<std::array<std::array<double,5>, dim> > large_array;
-
+  // Compute the values (and derivatives, if necessary) of all 1D polynomials
+  // at this evaluation point. We need to compute dim*n_polynomials
+  // evaluations, involving an evaluation of each polynomial for each
+  // coordinate direction. Once we have those values, we perform the
+  // multiplications for the tensor product in the arbitrary dimension.
   const unsigned int n_polynomials = polynomials.size();
-  if (n_polynomials > 20)
-    {
-      large_array.resize(n_polynomials);
-      v = &large_array[0];
-    }
+  boost::container::small_vector<std::array<std::array<double,5>,dim>, 20> values_1d(n_polynomials);
+  if (n_values_and_derivatives == 1)
+    for (unsigned int i=0; i<n_polynomials; ++i)
+      for (unsigned int d=0; d<dim; ++d)
+        values_1d[i][d][0] = polynomials[i].value(p(d));
   else
-    v = &small_array[0];
-
-  for (unsigned int i=0; i<n_polynomials; ++i)
-    for (unsigned int d=0; d<dim; ++d)
-      polynomials[i].value(p(d), n_values_and_derivatives, &v[i][d][0]);
-
-  for (unsigned int i=0; i<n_tensor_pols; ++i)
-    {
-      // first get the
-      // one-dimensional indices of
-      // this particular tensor
-      // product polynomial
-      unsigned int indices[dim];
-      compute_index (i, indices);
-
-      if (update_values)
-        {
-          values[i] = 1;
-          for (unsigned int x=0; x<dim; ++x)
-            values[i] *= v[indices[x]][x][0];
-        }
-
-      if (update_grads)
-        for (unsigned int d=0; d<dim; ++d)
+    for (unsigned int i=0; i<n_polynomials; ++i)
+      for (unsigned d=0; d<dim; ++d)
+        polynomials[i].value(p(d), n_values_and_derivatives, &values_1d[i][d][0]);
+
+  unsigned int indices[3];
+  unsigned int ind=0;
+  for (indices[2]=0; indices[2]<(dim>2?n_polynomials:1); ++indices[2])
+    for (indices[1]=0; indices[1]<(dim>1?n_polynomials:1); ++indices[1])
+      if (n_values_and_derivatives == 1)
+        for (indices[0]=0; indices[0]<n_polynomials; ++indices[0], ++ind)
           {
-            grads[i][d] = 1.;
-            for (unsigned int x=0; x<dim; ++x)
-              grads[i][d] *= v[indices[x]][x][d==x];
+            double value = values_1d[indices[0]][0][0];
+            for (unsigned int d=1; d<dim; ++d)
+              value *= values_1d[indices[d]][d][0];
+            values[index_map_inverse[ind]] = value;
           }
+      else
+        for (indices[0]=0; indices[0]<n_polynomials; ++indices[0], ++ind)
+          {
+            unsigned int i = index_map_inverse[ind];
 
-      if (update_grad_grads)
-        for (unsigned int d1=0; d1<dim; ++d1)
-          for (unsigned int d2=0; d2<dim; ++d2)
-            {
-              grad_grads[i][d1][d2] = 1.;
-              for (unsigned int x=0; x<dim; ++x)
-                {
-                  unsigned int derivative=0;
-                  if (d1==x) ++derivative;
-                  if (d2==x) ++derivative;
+            if (update_values)
+              {
+                double value = values_1d[indices[0]][0][0];
+                for (unsigned int x=1; x<dim; ++x)
+                  value *= values_1d[indices[x]][x][0];
+                values[i] = value;
+              }
 
-                  grad_grads[i][d1][d2]
-                  *= v[indices[x]][x][derivative];
+            if (update_grads)
+              for (unsigned int d=0; d<dim; ++d)
+                {
+                  double grad = 1.;
+                  for (unsigned int x=0; x<dim; ++x)
+                    grad *= values_1d[indices[x]][x][(d==x)?1:0];
+                  grads[i][d] = grad;
                 }
-            }
 
-      if (update_3rd_derivatives)
-        for (unsigned int d1=0; d1<dim; ++d1)
-          for (unsigned int d2=0; d2<dim; ++d2)
-            for (unsigned int d3=0; d3<dim; ++d3)
-              {
-                third_derivatives[i][d1][d2][d3] = 1.;
-                for (unsigned int x=0; x<dim; ++x)
+            if (update_grad_grads)
+              for (unsigned int d1=0; d1<dim; ++d1)
+                for (unsigned int d2=0; d2<dim; ++d2)
                   {
-                    unsigned int derivative=0;
-                    if (d1==x) ++derivative;
-                    if (d2==x) ++derivative;
-                    if (d3==x) ++derivative;
-
-                    third_derivatives[i][d1][d2][d3]
-                    *= v[indices[x]][x][derivative];
+                    double der2 = 1.;
+                    for (unsigned int x=0; x<dim; ++x)
+                      {
+                        unsigned int derivative=0;
+                        if (d1==x) ++derivative;
+                        if (d2==x) ++derivative;
+
+                        der2 *= values_1d[indices[x]][x][derivative];
+                      }
+                    grad_grads[i][d1][d2] = der2;
                   }
-              }
 
-      if (update_4th_derivatives)
-        for (unsigned int d1=0; d1<dim; ++d1)
-          for (unsigned int d2=0; d2<dim; ++d2)
-            for (unsigned int d3=0; d3<dim; ++d3)
-              for (unsigned int d4=0; d4<dim; ++d4)
-                {
-                  fourth_derivatives[i][d1][d2][d3][d4] = 1.;
-                  for (unsigned int x=0; x<dim; ++x)
+            if (update_3rd_derivatives)
+              for (unsigned int d1=0; d1<dim; ++d1)
+                for (unsigned int d2=0; d2<dim; ++d2)
+                  for (unsigned int d3=0; d3<dim; ++d3)
                     {
-                      unsigned int derivative=0;
-                      if (d1==x) ++derivative;
-                      if (d2==x) ++derivative;
-                      if (d3==x) ++derivative;
-                      if (d4==x) ++derivative;
-
-                      fourth_derivatives[i][d1][d2][d3][d4]
-                      *= v[indices[x]][x][derivative];
+                      double der3 = 1.;
+                      for (unsigned int x=0; x<dim; ++x)
+                        {
+                          unsigned int derivative=0;
+                          if (d1==x) ++derivative;
+                          if (d2==x) ++derivative;
+                          if (d3==x) ++derivative;
+
+                          der3 *= values_1d[indices[x]][x][derivative];
+                        }
+                      third_derivatives[i][d1][d2][d3] = der3;
                     }
-                }
-    }
+
+            if (update_4th_derivatives)
+              for (unsigned int d1=0; d1<dim; ++d1)
+                for (unsigned int d2=0; d2<dim; ++d2)
+                  for (unsigned int d3=0; d3<dim; ++d3)
+                    for (unsigned int d4=0; d4<dim; ++d4)
+                      {
+                        double der4 = 1.;
+                        for (unsigned int x=0; x<dim; ++x)
+                          {
+                            unsigned int derivative=0;
+                            if (d1==x) ++derivative;
+                            if (d2==x) ++derivative;
+                            if (d3==x) ++derivative;
+                            if (d4==x) ++derivative;
+
+                            der4 *= values_1d[indices[x]][x][derivative];
+                          }
+                        fourth_derivatives[i][d1][d2][d3][d4] = der4;
+                      }
+          }
 }