From: Martin Kronbichler 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=eb8597097dda2fb5da6084f6e1859b0c6e2c60ef;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 #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -245,6 +246,7 @@ compute (const Point &p, std::vector > &third_derivatives, std::vector > &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 &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 &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, dim> *v; - std::array, dim>, 20> small_array; - std::vector, 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,dim>, 20> values_1d(n_polynomials); + if (n_values_and_derivatives == 1) + for (unsigned int i=0; i2?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]