Using Eigen::VectorXd (Eigen 3.3.4) as a state typ

During my work, it would be a requirement for me to use Eigen::VectorXcd as state type, to solve a huge linear ODE system. In that project, the matrix in the ODE system is sparse. Multiplying a it with a dense vector can be computed in parallel in a simple way using Eigen. However, I have faced some problems in that situation that I will tell in details below. Currently I am applying a not-so efficient solution, namely, I use arma::cx_vec as state type, and declare the matrix corresponding the ode arma::cx_mat, which is a dense matrix type. Unfortunately, sparse-matrix-dense vector multiplication cannot be parallelized with Armadillo since it supports only the single-threaded version of SuperLU. According to my experience, in my specific case, sparse matrix-dense vector multiplication is three times faster than dense matrix-dense vector multiplication. Since a single running lasts for tens or hundreds of hours even when using 20 CPU cores simultaneously, it would be highly required for me to be able to use the Eigen library along with boost::odeint to gain a factor of three in runtime.

For testing purposes, I tried to use Eigen::VectorXd as a state type. I have also found very similar situations, see: Using Boost::odeint with Eigen::Matrix as a state vector or using several eigen matrices as statetypes in boost/odeint. However, none of the solutions worked for me. In my real project, I have to use the Bulirsch-Stoer method to integrate my ODE system since high precision and therefore error control is necessary for me. Below I will post my test code, corresponding to Eigen::VectorXd state type. As you can see, for first, I followed the steps described in the boost::odeint documentation: I defined my state vector resizeable and overloaded the mathematical operators used by the vector space algebra of the controlled stepper. After these steps, the code still did not compile, overloading of the addition operator for adding vectors to scalars was also necessary. Here is my code:

#include <iostream>
#include <vector>
#include <complex>

#include <Eigen/Eigen>
#include <boost/numeric/odeint.hpp>

namespace boost {
namespace numeric {
namespace odeint {

template<>
struct is_resizeable<Eigen::VectorXd>
{
    typedef boost::true_type type;
    static const bool value = type::value;
};

template<>
struct same_size_impl<Eigen::VectorXd, Eigen::VectorXd>
{
    static bool same_size(const Eigen::VectorXd& v1, const Eigen::VectorXd& v2)
    {
        std::cout << "same_size: " << v1.transpose() << " and " << v2.transpose() << std::endl;
        std::cout << "same_size: " << v1.rows() << " == " << v2.rows() << std::endl;
        return v1.rows() == v2.rows();
    }
};

template<>
struct resize_impl<Eigen::VectorXd, Eigen::VectorXd>
{
    static void resize(Eigen::VectorXd& v1, const Eigen::VectorXd& v2)
    {
        std::cout << "resize_impl: " << v1.transpose() << " and " << v2.transpose() << std::endl;
        std::cout << "resize_impl: (" << v1.rows() << " x " << v1.cols() << ") to (" << v2.rows() << " x " << v2.cols() << ")" << std::endl;
        v1.conservativeResize(v2.size());
        std::cout << "resize_impl: new size: " << v1.rows() << " x " << v1.cols() << std::endl;
    }
};

template<>
struct vector_space_norm_inf<Eigen::VectorXd>
{
    typedef double result_type;
    double operator()(const Eigen::VectorXd& v) const
    {
        std::cout << "vector_space_norm_inf(" << v.transpose() << ") = " << v.lpNorm<Eigen::Infinity>() << std::endl;
        return v.lpNorm<Eigen::Infinity>();
    }
};

}
}
}

namespace Eigen {
//VectorXd operator*(const double& a, const VectorXd& v)
//{
//  std::cout << a << "*(" << v.transpose() << ") = " << a*(v.transpose()) << std::endl;
//  return a*v;
//}
//
//VectorXd operator*(const VectorXd& v, const double& a)
//{
//  std::cout << "(" << v.transpose() << ")*" << v.transpose() << " = " << a*(v.transpose()) << std::endl;
//  return a*v;
//}
//
//VectorXd& operator*=(const double& a, VectorXd& v)
//{
//  std::cout << "operator*=(const double& a, VectorXd& v): (" << v.transpose() << ") *= " << a << ": ";
//  v *= a;
//  std::cout << v.transpose() << std::endl;
//  return v;
//}
//
//VectorXd& operator*=(VectorXd& v, const double& a)
//{
//  std::cout << "operator*=(VectorXd& v, const double& a): (" << v.transpose() << ") *= " << a << ": ";
//  v *= a;
//  std::cout << v.transpose() << std::endl;
//  return v;
//}
//
//VectorXd& operator+=(VectorXd& v1, const VectorXd& v2)
//{
//  std::cout << "(" << v1.transpose() << ") += (" << v2.transpose() << "): ";
//  v1 += v2;
//  std::cout << v1.transpose() << std::endl;
//  return v1;
//}

VectorXd operator/(const VectorXd& v1, const VectorXd& v2)
{
    const char* sizeeq;
    if(v1.size() == v2.size())
        sizeeq = "True";
    else
        sizeeq = "FALSE!";

    std::cout << "operator/: size check: " << v1.size() << " ?= " << v2.size() << ": " << sizeeq << std::endl;
    VectorXd result(v1.cwiseQuotient(v2));

    std::cout << "(" << v1.transpose() << ") / (" << v2.transpose() << ") = " << result.transpose() << std::endl;
    return result;
}

VectorXd operator+(double d, const VectorXd& v)
{
    std::cout << d << " + (" << v.transpose() << ") = " << d+v.transpose() << std::endl;
    return d+v;
}

VectorXd operator+(const VectorXd& v, double d)
{
    std::cout << "(" << v.transpose() << ") + " << d << " = " << d+v.transpose() << std::endl;
    return d+v;
}

VectorXd abs(const VectorXd& v)
{
    std::cout << "|(" << v.transpose() << ")| = " << v.cwiseAbs().transpose() << std::endl;
    return v.cwiseAbs();
}
}

typedef Eigen::VectorXd state_type;
typedef Eigen::MatrixXd coeff_matrix_type;

class lin_diff_eq
{
    public:
        lin_diff_eq(coeff_matrix_type gam): m_gam(gam) {}
        coeff_matrix_type get_matrix() const { return m_gam; }

        void operator()(const state_type& x, state_type& dxdt, const double t)
        {
            std::cout << "size: " << x.size() << std::endl;
            std::cout << x.transpose() << std::endl;
            dxdt = m_gam*x;
        }

    private:
        coeff_matrix_type m_gam;
};

using namespace boost::numeric::odeint;

int main()
{
    coeff_matrix_type A = coeff_matrix_type(2,2);
    state_type x = state_type(2);
    double epsabs = 1e-4;
    double epsrel = 1e-4;
    A << 0, 1,
        -2, 0;
    x << 1, 0;
    lin_diff_eq lde(A);
    std::cout << lde.get_matrix() << std::endl;

    integrate_adaptive(bulirsch_stoer<state_type, double, state_type, double, vector_space_algebra>(epsabs, epsrel), lin_diff_eq(A), x, 0.0, 1.0, 0.1);
}

All these steps made my code compile, however, during runtime my program crashes with segmentation fault. As you can see, I overloaded all the operators such that they produce some kind of a very simple log, I couldn't figure out how to produce a more detailed log to locate exactly where and why the segmentation fault occurs. According to the boost::odeint documentation, this can happen during runtime when the vector space algebra doesn't know how to resize the vectors, but I already defined how to do it. On my computer, it can be seen that boost::odeint generates some internal states, resizes them correctly, computes elementwise abs() correctly, but after that it crashes.

The next thing I tried was to use the support to Eigen provided by the boost::odeint developers. In this case, I did not use any of the template specializations and the operator overloads provided above, but instead, I included

boost/numeric/odeint/external/eigen/eigen_algebra.hpp

after including odeint.hpp. The code did not compile, I got the following error messages:

/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|35|error: ‘scalar_add_op’ in namespace ‘Eigen::internal’ does not name a type|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|35|error: expected template-argument before ‘<’ token|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|35|error: expected ‘>’ before ‘<’ token|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|37|error: wrong number of template arguments (1, should be 2)|
/usr/local/include/Eigen/src/Core/util/ForwardDeclarations.h|91|error: provided for ‘template<class UnaryOp, class MatrixType> class Eigen::CwiseUnaryOp’|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|38|error: expected ‘::’ before ‘operator’|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|38|error: expected identifier before ‘operator’|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp||In function ‘const int Eigen::operator+(const Eigen::MatrixBase<Derived>&, const typename Eigen::internal::traits<T>::Scalar&)’:|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|41|error: ‘scalar_add_op’ in namespace ‘Eigen::internal’ does not name a type|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|41|error: expected template-argument before ‘<’ token|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|41|error: expected ‘>’ before ‘<’ token|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|43|error: wrong number of template arguments (1, should be 2)|
/usr/local/include/Eigen/src/Core/util/ForwardDeclarations.h|91|error: provided for ‘template<class UnaryOp, class MatrixType> class Eigen::CwiseUnaryOp’|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|43|error: expected ‘::’ before ‘(’ token|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|43|error: expected identifier before ‘(’ token|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|43|error: ‘scalar_add_op’ is not a member of ‘Eigen::internal’|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|44|error: expected ‘(’ before ‘>’ token|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|50|error: ‘scalar_add_op’ in namespace ‘Eigen::internal’ does not name a type|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|50|error: expected template-argument before ‘<’ token|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|50|error: expected ‘>’ before ‘<’ token|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|52|error: wrong number of template arguments (1, should be 2)|
/usr/local/include/Eigen/src/Core/util/ForwardDeclarations.h|91|error: provided for ‘template<class UnaryOp, class MatrixType> class Eigen::CwiseUnaryOp’|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|53|error: expected ‘::’ before ‘operator’|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|53|error: expected identifier before ‘operator’|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp||In function ‘const int Eigen::operator+(const typename Eigen::internal::traits<T>::Scalar&, const Eigen::MatrixBase<Derived>&)’:|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|56|error: ‘scalar_add_op’ in namespace ‘Eigen::internal’ does not name a type|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|56|error: expected template-argument before ‘<’ token|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|56|error: expected ‘>’ before ‘<’ token|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|58|error: wrong number of template arguments (1, should be 2)|
/usr/local/include/Eigen/src/Core/util/ForwardDeclarations.h|91|error: provided for ‘template<class UnaryOp, class MatrixType> class Eigen::CwiseUnaryOp’|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|58|error: expected ‘::’ before ‘(’ token|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|58|error: expected identifier before ‘(’ token|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|58|error: ‘scalar_add_op’ is not a member of ‘Eigen::internal’|
/usr/local/include/boost/numeric/odeint/external/eigen/eigen_algebra.hpp|59|error: expected ‘(’ before ‘>’ token|
||=== Build failed: 32 error(s), 0 warning(s) (0 minute(s), 7 second(s)) ===|

I also tried to include

boost/numeric/odeint/external/eigen/eigen.hpp

instead of eigen_algebra.hpp and received the same error messages when including eigen_algebra.hpp. It seems like Boost::odeint v1.65.1 is not compatible with Eigen v3.3.4. My question is if there is an (un)official header file that provides a support of Eigen 3.3.4 to Boost::odeint 1.65.1, or did I make something wrong during template specializations and/or operator overloading for the vector space algebra?

I also found a solution for "typecasting" std std::vector<> to Eigen::VectorXd and vice versa, but I think that using std::vector as a state stype, and then internally "casting it" to Eigen::VectorXd, then "casting the result back" to std::vector<> would result in a lower performance than that of mentioned in the introduction since my vectors are huge and ten tousands or more accepted steps are required to perform an integration. I hope my question was accurate enough and the community can kindly help me solving my problem.