I'm writing a numerical methods for a specific ODE of the form
y'(t) = A*y(t) + g(t,y) in C++, where A is a matrix and g is a vector function, i.e. a function that takes as input the time and a vector and return a vector. I'm using the library Eigen.
I'm writing my own class, but I can't understand how to pass the vector function g(t,y) as an argument in the class definition, which is
integrator(double fin_time, int time_steps, MatrixXd rhs, VectorXd g, VectorXd in_data)
How can I do that with Eigen?
When the function is scalar, I usually do double (*f)(double t, double y), but here is different of course.
Instead of using function pointers I prefer to use proper C++11 function objects. The function g(t, y) would then have the following type;
std::function<VectorXd (double, const VectorXd&)> g;
This type could be used in your class constructor.
As you mentioned, the typical size of your vectors is n=8000. For that size you can either use dynamic or fixed size Eigen vectors. If you use dynamic vectors you should avoid dynamic resizing. This can be achieved by allocating the correct size when defining the variable:
VectorXd vec(n);
Since you know the size n at compile time, you could also use a template parameter to fix the size:
template<int n>
class SpecialODE {
using Vector = Eigen::Matrix<double, n, 1>
};
Loop unrolling and SIMD vectorization is done by Eigen and the compiler in both cases. I would not expect to see any major differences in performance for larger matrix sizes.