Armadillo matrix as the state type

[patryk-kubiczek]

I am trying to use Armadillo complex matrix as the state type with a dense output stepper. Therefore, I am using vector_space_algebra with an instantiated vector_space_norm_inf. I put the following code in the header

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

using matrix_t = arma::cx_mat;

namespace boost {
    namespace numeric {
        namespace odeint {
            template<>
            struct vector_space_norm_inf<matrix_t> {
                typedef double return_type;
                return_type operator()(matrix_t &u) const {
                    return arma::norm(u, "inf");
                }
            };
        }}}

Then in the code I define the stepper type as

using stepper_t = dense_output_runge_kutta<controlled_runge_kutta<runge_kutta_dopri5<matrix_t, double, matrix_t, double, vector_space_algebra>>>;

However, during the compilation I get the following error

error: no matching function for call to ‘boost::numeric::odeint::vector_space_algebra::norm_inf(arma::Mat<std::complex<double> >&) return algebra.norm_inf( x_err );

I don't understand this, since I believe norm_inf has been specialized for Armadillo complex matrix. The problem occurs also for real matrices.

I would appreciate your help.

[mariomulansky]

For some reason, we called the type result_type, not return_type. Try using this name in the typedef and see if it works. Otherwise, your approach is correct.

[patryk-kubiczek]

Thank you for the hint, but it did not help. Actually, I am still having a similar error message

error: no match for call to ‘(boost::numeric::odeint::vector_space_norm_inf<arma::Mat<std::complex<double> > >) (const arma::Mat<std::complex<double> >&)’
         return n( s );

followed by

note: candidate: boost::numeric::odeint::vector_space_norm_inf<arma::Mat<std::complex<double> > >::result_type boost::numeric::odeint::vector_space_norm_inf<arma::Mat<std::complex<double> > >::operator()(matrix_t&) const <near match>
         result_type operator()(matrix_t &u) const {

note:   conversion of argument 1 would be ill-formed:

error: binding ‘const arma::Mat<std::complex<double> >’ to reference of type ‘matrix_t& {aka arma::Mat<std::complex<double> >&}’ discards qualifiers
         return n( s );

Could it be that the argument of operator() should be a (const) value and not reference? When the argument is matrix_t u there is no deduction error but rather Armadillo error

error: static assertion failed: error: incorrect or unsupported type
   arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));

which probably doesn't tell you much unfortunately.

[mariomulansky]

Did you try defining it as a const ref? This should be the right way.

[patryk-kubiczek]

Indeed, const ref seems to be the correct argument of operator(), i.e. now I have return_type operator()(const matrix_t &u). However, then some assertion in Armadillo fails and the code still does not compile:

/usr/include/armadillo_bits/Mat_meat.hpp: In instantiation of ‘const arma::Mat<eT>& arma::Mat<eT>::operator=(const arma::eGlue<T1, T2, eglue_type>&) [with T1 = arma::mtOp<double, arma::Mat<std::complex<double> >, arma::op_abs>; T2 = arma::eOp<arma::eOp<arma::eGlue<arma::eOp<arma::mtOp<double, arma::Mat<std::complex<double> >, arma::op_abs>, arma::eop_scalar_times>, arma::eOp<arma::mtOp<double, arma::Mat<std::complex<double> >, arma::op_abs>, arma::eop_scalar_times>, arma::eglue_plus>, arma::eop_scalar_times>, arma::eop_scalar_plus>; eglue_type = arma::eglue_div; eT = std::complex<double>]’:
(...)
/usr/include/armadillo_bits/Mat_meat.hpp:5117:3: error: static assertion failed: error: incorrect or unsupported type
   arma_type_check(( is_same_type< eT, typename T1::elem_type >::no ));
/usr/include/armadillo_bits/Mat_meat.hpp:5118:3: error: static assertion failed: error: incorrect or unsupported type
   arma_type_check(( is_same_type< eT, typename T2::elem_type >::no ));

I am using gcc version 5.4.0 and Armadillo version 6.500.