Re: [eigen] calculating on ranges

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

Hi,

Yes, Eigen's lazy evaluation means that if you implement this with
Eigen it will also evaluate the function only on the range that's
needed.

Compared to your C function, I don't know if you would prefer the
Eigen way or not. Normally, Eigen would bring the advantage of
vectorization, except that here you use exp and that is not yet
vectorized (it could be if that were important. we would use Pade
approximants if there are no SSE instructions for that.). The other
advantage of using Eigen is that it will automatically unroll the loop
if the size (here 10) is known at compile-time and the total cost
isn't too big (you can tune EIGEN_UNROLLING_LIMIT).

Here is a basic Eigen implementation:

#include <Eigen/Array>

template<typename VectorType>
void addGaussianToCurve(
const MatrixBase<VectorType>& x,
const MatrixBase<VectorType>& y,
typename VectorType::Scalar a,
typename VectorType::Scalar b,
typename VectorType::Scalar c)
{
  y += a * ( -(x.cwise() - b).cwise().square() / c ).cwise().exp();
}

Now if you have vectors VectorXd x,y; you can call your function like this:
addGaussianToCurve(
x.segment(start,length),
y.segment(start,length),
a,b,c);

Now if you know length at compile time you can do:
addGaussianToCurve(
x.segment<length>(start),
y.segment<length>(start),
a,b,c);
And it will automatically unroll the loop if appropriate (see above).

Again if vectorization matters it's no big deal, we just need a packet
version of exp(). You'd then benefit from it automatically.

Cheers,
Benoit