From ea0cc1d135d3023e5cbce60637e950a14572e77e Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sun, 11 Oct 2020 22:35:22 +0200 Subject: [PATCH] Make a polynomial function inline to allow templated use --- include/deal.II/base/polynomial.h | 161 ++++++++++++++++++++++++++++-- source/base/polynomial.cc | 142 -------------------------- 2 files changed, 153 insertions(+), 150 deletions(-) diff --git a/include/deal.II/base/polynomial.h b/include/deal.II/base/polynomial.h index 1ca3d93379..1c708c8cd5 100644 --- a/include/deal.II/base/polynomial.h +++ b/include/deal.II/base/polynomial.h @@ -93,9 +93,11 @@ namespace Polynomials /** * Return the value of this polynomial at the given point. * - * This function uses the Horner scheme for numerical stability of the - * evaluation for polynomials in the coefficient form or the product of - * terms involving the roots if that representation is used. + * This function uses the most numerically stable evaluation + * algorithm for the provided form of the polynomial. If the + * polynomial is in the product form of roots, the evaluation is + * based on products of the form (x - x_i), whereas the Horner + * scheme is used for polynomials in the coefficient form. */ number value(const number x) const; @@ -120,14 +122,22 @@ namespace Polynomials * determined by @p n_derivatives and @p values has to provide sufficient * space for @p n_derivatives + 1 values. * - * This function uses the Horner scheme for numerical stability of the - * evaluation for polynomials in the coefficient form or the product of - * terms involving the roots if that representation is used. + * This function uses the most numerically stable evaluation + * algorithm for the provided form of the polynomial. If the + * polynomial is in the product form of roots, the evaluation is + * based on products of the form (x - x_i), whereas the Horner + * scheme is used for polynomials in the coefficient form. + * + * The template type `Number2` must implement arithmetic + * operations such as additions or multiplication with the type + * `number` of the polynomial, and must be convertible from + * `number` by `operator=`. */ + template void - value(const number x, + value(const Number2 x, const unsigned int n_derivatives, - number * values) const; + Number2 * values) const; /** * Degree of the polynomial. This is the degree reflected by the number of @@ -807,6 +817,141 @@ namespace Polynomials + template + template + inline void + Polynomial::value(const Number2 x, + const unsigned int n_derivatives, + Number2 * values) const + { + // evaluate Lagrange polynomial and derivatives + if (in_lagrange_product_form == true) + { + // to compute the value and all derivatives of a polynomial of the + // form (x-x_1)*(x-x_2)*...*(x-x_n), expand the derivatives like + // automatic differentiation does. + const unsigned int n_supp = lagrange_support_points.size(); + const number weight = lagrange_weight; + switch (n_derivatives) + { + default: + values[0] = 1.; + for (unsigned int d = 1; d <= n_derivatives; ++d) + values[d] = 0.; + for (unsigned int i = 0; i < n_supp; ++i) + { + const Number2 v = x - lagrange_support_points[i]; + + // multiply by (x-x_i) and compute action on all derivatives, + // too (inspired from automatic differentiation: implement the + // product rule for the old value and the new variable 'v', + // i.e., expand value v and derivative one). since we reuse a + // value from the next lower derivative from the steps before, + // need to start from the highest derivative + for (unsigned int k = n_derivatives; k > 0; --k) + values[k] = (values[k] * v + values[k - 1]); + values[0] *= v; + } + // finally, multiply by the weight in the Lagrange + // denominator. Could be done instead of setting values[0] = 1 + // above, but that gives different accumulation of round-off + // errors (multiplication is not associative) compared to when we + // computed the weight, and hence a basis function might not be + // exactly one at the center point, which is nice to have. We also + // multiply derivatives by k! to transform the product p_n = + // p^(n)(x)/k! into the actual form of the derivative + { + number k_factorial = 1; + for (unsigned int k = 0; k <= n_derivatives; ++k) + { + values[k] *= k_factorial * weight; + k_factorial *= static_cast(k + 1); + } + } + break; + + // manually implement case 0 (values only), case 1 (value + first + // derivative), and case 2 (up to second derivative) since they + // might be called often. then, we can unroll the inner loop and + // keep the temporary results as local variables to help the + // compiler with the pointer aliasing analysis. + case 0: + { + Number2 value = 1.; + for (unsigned int i = 0; i < n_supp; ++i) + { + const Number2 v = x - lagrange_support_points[i]; + value *= v; + } + values[0] = weight * value; + break; + } + + case 1: + { + Number2 value = 1.; + Number2 derivative = 0.; + for (unsigned int i = 0; i < n_supp; ++i) + { + const Number2 v = x - lagrange_support_points[i]; + derivative = derivative * v + value; + value *= v; + } + values[0] = weight * value; + values[1] = weight * derivative; + break; + } + + case 2: + { + Number2 value = 1.; + Number2 derivative = 0.; + Number2 second = 0.; + for (unsigned int i = 0; i < n_supp; ++i) + { + const Number2 v = x - lagrange_support_points[i]; + second = second * v + derivative; + derivative = derivative * v + value; + value *= v; + } + values[0] = weight * value; + values[1] = weight * derivative; + values[2] = static_cast(2) * weight * second; + break; + } + } + return; + } + + Assert(coefficients.size() > 0, ExcEmptyObject()); + + // if derivatives are needed, then do it properly by the full + // Horner scheme + const unsigned int m = coefficients.size(); + std::vector a(coefficients.size()); + std::copy(coefficients.begin(), coefficients.end(), a.begin()); + unsigned int j_factorial = 1; + + // loop over all requested derivatives. note that derivatives @p{j>m} are + // necessarily zero, as they differentiate the polynomial more often than + // the highest power is + const unsigned int min_valuessize_m = std::min(n_derivatives + 1, m); + for (unsigned int j = 0; j < min_valuessize_m; ++j) + { + for (int k = m - 2; k >= static_cast(j); --k) + a[k] += x * a[k + 1]; + values[j] = static_cast(j_factorial) * a[j]; + + j_factorial *= j + 1; + } + + // fill higher derivatives by zero + for (unsigned int j = min_valuessize_m; j <= n_derivatives; ++j) + values[j] = 0.; + } + + + template template inline void diff --git a/source/base/polynomial.cc b/source/base/polynomial.cc index 7bf897c29b..9f9efbeb8d 100644 --- a/source/base/polynomial.cc +++ b/source/base/polynomial.cc @@ -106,148 +106,6 @@ namespace Polynomials - template - void - Polynomial::value(const number x, - const unsigned int n_derivatives, - number * values) const - { - // evaluate Lagrange polynomial and derivatives - if (in_lagrange_product_form == true) - { - // to compute the value and all derivatives of a polynomial of the - // form (x-x_1)*(x-x_2)*...*(x-x_n), expand the derivatives like - // automatic differentiation does. - const unsigned int n_supp = lagrange_support_points.size(); - const number weight = lagrange_weight; - switch (n_derivatives) - { - default: - values[0] = 1; - for (unsigned int d = 1; d <= n_derivatives; ++d) - values[d] = 0; - for (unsigned int i = 0; i < n_supp; ++i) - { - const number v = x - lagrange_support_points[i]; - - // multiply by (x-x_i) and compute action on all derivatives, - // too (inspired from automatic differentiation: implement the - // product rule for the old value and the new variable 'v', - // i.e., expand value v and derivative one). since we reuse a - // value from the next lower derivative from the steps before, - // need to start from the highest derivative - for (unsigned int k = n_derivatives; k > 0; --k) - values[k] = (values[k] * v + values[k - 1]); - values[0] *= v; - } - // finally, multiply by the weight in the Lagrange - // denominator. Could be done instead of setting values[0] = 1 - // above, but that gives different accumulation of round-off - // errors (multiplication is not associative) compared to when we - // computed the weight, and hence a basis function might not be - // exactly one at the center point, which is nice to have. We also - // multiply derivatives by k! to transform the product p_n = - // p^(n)(x)/k! into the actual form of the derivative - { - number k_factorial = 1; - for (unsigned int k = 0; k <= n_derivatives; ++k) - { - values[k] *= k_factorial * weight; - k_factorial *= static_cast(k + 1); - } - } - break; - - // manually implement case 0 (values only), case 1 (value + first - // derivative), and case 2 (up to second derivative) since they - // might be called often. then, we can unroll the inner loop and - // keep the temporary results as local variables to help the - // compiler with the pointer aliasing analysis. - case 0: - { - number value = 1; - for (unsigned int i = 0; i < n_supp; ++i) - { - const number v = x - lagrange_support_points[i]; - value *= v; - } - values[0] = weight * value; - break; - } - - case 1: - { - number value = 1; - number derivative = 0; - for (unsigned int i = 0; i < n_supp; ++i) - { - const number v = x - lagrange_support_points[i]; - derivative = derivative * v + value; - value *= v; - } - values[0] = weight * value; - values[1] = weight * derivative; - break; - } - - case 2: - { - number value = 1; - number derivative = 0; - number second = 0; - for (unsigned int i = 0; i < n_supp; ++i) - { - const number v = x - lagrange_support_points[i]; - second = second * v + derivative; - derivative = derivative * v + value; - value *= v; - } - values[0] = weight * value; - values[1] = weight * derivative; - values[2] = static_cast(2) * weight * second; - break; - } - } - return; - } - - Assert(coefficients.size() > 0, ExcEmptyObject()); - - // if we only need the value, then call the other function since that is - // significantly faster (there is no need to allocate and free memory, - // which is really expensive compared to all the other operations!) - if (n_derivatives == 0) - { - values[0] = value(x); - return; - } - - // if there are derivatives needed, then do it properly by the full Horner - // scheme - const unsigned int m = coefficients.size(); - std::vector a(coefficients); - unsigned int j_factorial = 1; - - // loop over all requested derivatives. note that derivatives @p{j>m} are - // necessarily zero, as they differentiate the polynomial more often than - // the highest power is - const unsigned int min_valuessize_m = std::min(n_derivatives + 1, m); - for (unsigned int j = 0; j < min_valuessize_m; ++j) - { - for (int k = m - 2; k >= static_cast(j); --k) - a[k] += x * a[k + 1]; - values[j] = static_cast(j_factorial) * a[j]; - - j_factorial *= j + 1; - } - - // fill higher derivatives by zero - for (unsigned int j = min_valuessize_m; j <= n_derivatives; ++j) - values[j] = 0; - } - - - template void Polynomial::transform_into_standard_form() -- 2.39.5