I have the following function:
double
neville (double xx, size_t n, const double *x, const double *y, double *work);
which performs Lagrange interpolation at xx
using the n
points stored in x
and y
. The work
array has size 2 * n
. Since this is polynomial interpolation, n
is in the ballpark of ~5, very rarely more than 10.
This function is aggressively optimized, and supposed to be called in tight loops. Profiling suggests that heap allocating the work array in the loop is bad. Unfortunately, I'm supposed to package this into a function-like class, and clients must be unaware of the work array.
For now, I use a template integer argument for the degree and std::array
to avoid dynamic allocation of the work
array:
template <size_t n>
struct interpolator
{
double operator() (double xx) const
{
std::array<double, 2 * n> work;
size_t i = locate (xx); // not shown here, no performance impact
// due to clever tricks + nice calling patterns
return neville (xx, n, x + i, y + i, work.data ());
}
const double *x, *y;
};
It would have been possible to store the work array as a mutable member of the class, but operator()
is supposed to be used concurrently by several threads. This version is OK provided you know n
at compile time.
Now, I need the n
parameter to be specified at run time. I am wondering about something like this:
double operator() (double xx) const
{
auto work = static_cast<double*> (alloca (n * sizeof (double)));
...
Some bells ring when using alloca
: I'm of course going to have a cap on n
to avoid the alloca
call to overflow (anyway it's quite stupid to use degree 100 polynomial interpolation).
I'm quite unconfortable with the approach however:
- Am I missing some obvious danger of
alloca
? - Is there a better way to avoid heap allocation here ?