Improve basis function implementation #1338
Conversation
davidscn
added
the
enhancement
src/mapping/impl/BasisFunctions.hpp
Outdated
| using std::log; | ||
| using std::pow; | ||
| return 1.0 - 30.0 * pow(p, 2.0) - 10.0 * pow(p, 3.0) + 45.0 * pow(p, 4.0) - 6.0 * pow(p, 5.0) - 60.0 * log(pow(p, pow(p, 3.0))); | ||
| return 1.0 - 30.0 * math::pow_int<2>(p) - 10.0 * math::pow_int<3>(p) + 45.0 * math::pow_int<4>(p) - 6.0 * math::pow_int<5>(p) - math::pow_int<3>(p) * 60.0 * std::log(std::max(p, 1e-15)); |
There was a problem hiding this comment.
Instead of going for this overly complicated formula in the end log(p^(p^3)) I selected now a fixed cutoff for the log (1e-14). The overall product will be zero anyway, as well multiply it by (1e-14)^3.
|
It is obscure, but there actually is an overload that matches a double base and a integral exponent. Its overload 7 of the docs. So, the exponent should be written 2 instead of 2.0. |
I was indeed a bit puzzled about this version. The docs say "If any argument has integral type, it is cast to double." So, I guess it is still double? Anyway, the recursion version is still faster than the |
|
Ah I see. You lifted the exponent into the type system. Makes sense that it is faster. You did test using a Release build, right? |
Yes (although I would have expected a somewhat similar result for small exponents from the STL).
Yes. |
* Implement basis functions via recursion * Add unit test for pow method * Use math::NUMERICAL_ZERO_DIFFERENCE for tolerance * Apply suggestions from code review and add changelog entry
Main changes of this PR
Switches the basis function implementation from using
std::powto a pow variant relying on recursion.Motivation and additional information
The main problem is that
std::powoperates on adoubleexponent (there is nointvariant), but we use it in the basis function implementation only for integer exponents. The basis function evaluation routine is heavily performance critical. With the changes here, we get a considerable performance gain in both RBF implementations (Eigen and PETSc). Here some timings for the (expensive) computation of matrix A as well as the overall mapping computation:result-asmA.pdf
result-computet.pdf
Example: For the c6 Wendland RBF, we see a speed-up of ~4.
Also related to #1320
Author's checklist
make changelogif there are user-observable changes since the last release.make formatto ensure everything is formatted correctly.Reviewers' checklist