function pointer with Eigen

Question

I am pretty competent with Python, but I'm pretty new to C++ and things like pointers. I try to write some codes for solving ODE with the Eigen package for linear algebra (I will need to deal with lots of matrices later, so I plan to start with it). I have the following code for RK4 and they work:

#include "../eigen-eigen-b3f3d4950030/Eigen/Dense"
using namespace Eigen;

VectorXd Func(const VectorXd& a)
{ // equations for solving simple harmonic oscillator
    Vector2d ans;
    ans(0) = a(1);    //  dy/dt
    ans(1) = -a(0);   //  d2y/dt2
    return ans;
}

MatrixXd RK4(VectorXd Func(const VectorXd& y), const Ref<const VectorXd>& y0, double h, int step_num)
{
    MatrixXd y(step_num, y0.rows());
    y.row(0) = y0;

    for (int i=1; i<step_num; i++){
        VectorXd y_old = y.row(i-1).transpose();
        VectorXd k1 = h*Func(y_old);
        VectorXd k2 = h*Func(y_old+k1/2);
        VectorXd k3 = h*Func(y_old+k2/2);
        VectorXd k4 = h*Func(y_old+k3);
        VectorXd dy = (k1 + 2*k2 + 2*k3 + k4)/6;
        y.row(i) = y.row(i-1) + dy.transpose();
    }
    return y;
}

int main()
{
    Vector2d v1;
    v1(0) = 1.4;    v1(1) = -0.1;
    double h = 0.1;
    int step_num = 50;

    MatrixXd sol = RK4(Func,v1,h,step_num);
    return 0;
}

I have the following questions:

What's the meaning of & in the function argument? Pass by reference? I just copied the code from the official documentation, but I'm not too sure if I understand every bit in the function arguments of RK4 such as VectorXd Func(const VectorXd& y). Are there alternative ways of accepting Eigen::MatrixXd and functions which accept Eigen::MatrixXd as function arguments?
From what I understand, we cannot return a whole 2D array from a function, and what we are returning is just the first element of the array (correct me if I'm wrong). What about the Eigen::MatrixX? What are we actually passing/returning? The first element of the matrix, or a completely new object defined by the Eigen library?
I'm not sure if these codes are written efficiently. Are there anything I can do to optimize this part? (Just wondering if I have done anything that may significantly slow down the speed).

Thanks

a) please one question per question b) take a look [here](https://stackoverflow.com/questions/388242/the-definitive-c-book-guide-and-list) — 463035818_is_not_an_ai, Aug 14 '18 at 11:41
Regarding your first question: Yes you are right. A function declared as `ReturnType FunctionName(ArgType& arg);` Takes a reference to a variable of type `ArgType` called `arg`. Note: Because (specially) for bigger Objects passing the Argument by reference instead of copying them leads to better performance many functions use `const` References for their arguments. See: https://stackoverflow.com/questions/3694630/c-const-reference-before-vs-after-type-specifier — maperz, Aug 14 '18 at 12:24
Regarding performance, use compile-time vector/matrix types whenever possible for tiny objects (say <10), like `Vector2d`. — ggael, Aug 14 '18 at 14:01

score 1 · Accepted Answer · edited Jun 20 '20 at 09:12

Yes, & is pass-by-reference; The latter one is syntax for passing a function, that takes a vector by reference and returns a vector. An Eigen::Matrix should always be passed by reference. There are tons of ways to pass one function to another, the most idiomatic ones in C++ are probably template arguments and std::function.
You can't have multiple return arguments, but you can return a pair or a tuple or a Matrix object. RK4 returns a whole matrix.
The code is fairly efficient. If it was really performance-critical there might be a few things that could be optimized, but I would not worry for now.

The biggest point is that RK4 is very general and works with dynamically sized types, which are a lot more expensive than their statically sized counter parts (VectorXf vs Vector2d). But this would require you to create a specialized version for all dimensions you are interested in or to get the compiler to do it for you by using templates.

Generally: Read a good book to get you started.

function pointer with Eigen

1 Answers1

Linked