I came around some weird behavior concerning the eigen library and templated functions.
Maybe someone can explain to me, why the first version is not working, while the other 3 do. My guess would be the first case freeing up some local variable, but hopefully someone can enlighten me. Thanks in advance.
Here is the code:
Compiler-Explorer: https://compiler-explorer.com/z/r45xzE417
#include <concepts>
#include <iostream>
#include <Eigen/Core>
auto RungeKutta1_auto(const auto& f, const std::floating_point auto& h, const auto& y_n)
{
auto ret = y_n + h * f(y_n);
std::cout << ret.transpose() << std::endl;
return ret;
}
template<typename _Scalar, int _Rows, int _Cols>
auto RungeKutta1_template(const auto& f, const std::floating_point auto& h, const Eigen::Matrix<_Scalar, _Rows, _Cols>& y_n)
{
Eigen::Matrix<_Scalar, _Rows, _Cols> ret = y_n + h * f(y_n);
std::cout << ret.transpose() << std::endl;
return ret;
}
int main()
{
auto f = [](const Eigen::Matrix<double, 10, 1>& y) {
Eigen::Matrix<double, 10, 1> y_dot = 2 * y;
return y_dot;
};
auto g = [](const Eigen::Matrix<double, 10, 1>& y) {
return 2 * y;
};
std::cout << "RungeKutta1_auto(f, 0.05, y):" << std::endl;
Eigen::Matrix<double, 10, 1> y;
y << 0, 1, 2, 3, 4, 5, 6, 7, 8, 9;
y = RungeKutta1_auto(f, 0.05, y);
std::cout << y.transpose() << std::endl;
// Output
// 0 1.1 2.2 3.3 4.4 5.5 6.6 7.7 8.8 9.9
// 3.47627e-311 1 2 3 4 5 6.6 7 8 9
std::cout << "RungeKutta1_template(f, 0.05, y):" << std::endl;
y << 0, 1, 2, 3, 4, 5, 6, 7, 8, 9;
y = RungeKutta1_template(f, 0.05, y);
std::cout << y.transpose() << std::endl;
// Output
// 0 1.1 2.2 3.3 4.4 5.5 6.6 7.7 8.8 9.9
// 0 1.1 2.2 3.3 4.4 5.5 6.6 7.7 8.8 9.9
std::cout << "RungeKutta1_auto(g, 0.05, y):" << std::endl;
y << 0, 1, 2, 3, 4, 5, 6, 7, 8, 9;
y = RungeKutta1_auto(g, 0.05, y);
std::cout << y.transpose() << std::endl;
// Output
// 0 1.1 2.2 3.3 4.4 5.5 6.6 7.7 8.8 9.9
// 0 1.1 2.2 3.3 4.4 5.5 6.6 7.7 8.8 9.9
std::cout << "RungeKutta1_template(g, 0.05, y):" << std::endl;
y << 0, 1, 2, 3, 4, 5, 6, 7, 8, 9;
y = RungeKutta1_template(g, 0.05, y);
std::cout << y.transpose() << std::endl;
// Output
// 0 1.1 2.2 3.3 4.4 5.5 6.6 7.7 8.8 9.9
// 0 1.1 2.2 3.3 4.4 5.5 6.6 7.7 8.8 9.9
}
In the first version,
auto ret = y_n + h * f(y_n);
due to Eigen's expression templates this gives you an intermediate expression type, as opposed to a Matrix type. You would need to explicitly invoke eval() on it to force the lazy eval to occur (and in such converting the expression to a resulting Matrix type object); e.g.:
auto ret = (y_n + h * f(y_n)).eval();
In the other versions you are explicitly typing out the type of ret to be a Matrix type, meaning you will not end up storing intermediate expression template types.