Multivariate Function with Known Gradients - RStan

The solution to this issue is quite simple as I got know in the thread Barriers to New Developers. Basically, one should not expose a vector-valued function based on std::vector, but always based on Eigen::Matrix. Additionally, there always has to be an overloaded function definition for both double and var.

namespace stan {
    namespace math {

        inline double
        vecToDoubleTwo(const Eigen::Matrix<double, Eigen::Dynamic, 1>& m) {

            double val = 0;
            for (int i=0; i<m.size(); ++i) val += m(i) * m(i);

            return val;
        }


        inline var
        vecToDoubleTwo(const Eigen::Matrix<var, Eigen::Dynamic, 1>& m) {
            using Eigen::Matrix;

            Matrix<double, Eigen::Dynamic, 1> m_d(m.rows());
            for (int i = 0; i < m.size(); ++i)
                m_d(i) = m(i).val();

            double val = stan::math::vecToDouble(m_d);

            vari** varis
            = ChainableStack::memalloc_.alloc_array<vari*>(m.size());
            for (int i = 0; i < m.size(); ++i)
                varis[i] = m(i).vi_;

            double* gradients
            = ChainableStack::memalloc_.alloc_array<double>(m.size());
            for (int i = 0; i < m.size(); ++i)
                gradients[i] = 9.9 * m_d(i);

            return var(new precomputed_gradients_vari(val, m.size(),
                                                      varis, gradients));
        }

    }
}