I have the following code:

Code: Select all

        template <typename T>
        T my_dot (const Eigen::ArrayXd& a, const Eigen::ArrayX<T>& b) {
            const Eigen::array<const Eigen::IndexPair<int>, 1> product_dims = {Eigen::IndexPair<int>(0,0)};
            const Eigen::Tensor<double, 1> a_tensor = Eigen::TensorMap<const Eigen::Tensor<const double, 1>> (a.data(), a.size());;
            const Eigen::Tensor<T, 1> b_tensor = Eigen::TensorMap<const Eigen::Tensor<const T, 1>> (b.data(), b.size());
            const Eigen::Tensor<T, 0> ret = a_tensor.contract(b_tensor, product_dims);
            return ret();
        }

        template <typename T>
        T my_dot2 (const Eigen::ArrayXd& a, const Eigen::ArrayX<T>& b) {
            return (a*b).sum();
        }

my_dot is templated so as to allow automatic differentiation via stan. So the substitution I'm most interested in is T=stan::math::var, with is essentially a pair of doubles with funky arithmetic. (My understanding of how backwards autodiff is implemented is fairly thin, so this might be wrong).

My issue is that my_dot2 works perfectly with autodiff types, but my_dot won't compile.

Code: Select all

/usr/include/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:1555:15: error: no matching function for call to ‘Eigen::internal::gebp_traits<stan::math::var, double, false, false>::madd(Eigen::internal::gebp_kernel<double, stan::math::var, long int, Eigen::internal::blas_data_mapper<stan::math::var, long int, 0, 0>, 2, 4, false, false>::SLhsPacket&, Eigen::internal::gebp_kernel<double, stan::math::var, long int, Eigen::internal::blas_data_mapper<stan::math::var, long int, 0, 0>, 2, 4, false, false>::SRhsPacket&, Eigen::internal::gebp_kernel<double, stan::math::var, long int, Eigen::internal::blas_data_mapper<stan::math::var, long int, 0, 0>, 2, 4, false, false>::SAccPacket&, Eigen::internal::gebp_kernel<double, stan::math::var, long int, Eigen::internal::blas_data_mapper<stan::math::var, long int, 0, 0>, 2, 4, false, false>::SRhsPacket&)’
               straits.madd(A0,B_0,C0,B_0);
               ^~~~~~~
/usr/include/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:435:28: note: candidate: template<class LhsPacketType, class RhsPacketType, class AccPacketType> void Eigen::internal::gebp_traits<_LhsScalar, _RhsScalar, _ConjLhs, _ConjRhs>::madd(const LhsPacketType&, const RhsPacketType&, AccPacketType&, AccPacketType&) const [with LhsPacketType = LhsPacketType; RhsPacketType = RhsPacketType; AccPacketType = AccPacketType; _LhsScalar = stan::math::var; _RhsScalar = double; bool _ConjLhs = false; bool _ConjRhs = false]
   EIGEN_STRONG_INLINE void madd(const LhsPacketType& a, const RhsPacketType& b, AccPacketType& c, AccPacketType& tmp) const
                            ^~~~
/usr/include/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:435:28: note:   template argument deduction/substitution failed:
/usr/include/eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h:1555:15: note:   deduced conflicting types for parameter ‘AccPacketType’ (‘stan::math::var’ and ‘Eigen::internal::gebp_kernel<double, stan::math::var, long int, Eigen::internal::blas_data_mapper<stan::math::var, long int, 0, 0>, 2, 4, false, false>::SRhsPacket {aka double}’)
               straits.madd(A0,B_0,C0,B_0);
               ^~~~~~~

Or if I change my_dot by swapping the order of arguments:

Code: Select all

        template <typename T>
        T my_dot (const Eigen::ArrayXd& a, const Eigen::ArrayX<T>& b) {
            const Eigen::array<const Eigen::IndexPair<int>, 1> product_dims = {Eigen::IndexPair<int>(0,0)};
            const Eigen::Tensor<double, 1> a_tensor = Eigen::TensorMap<const Eigen::Tensor<const double, 1>> (a.data(), a.size());;
            const Eigen::Tensor<T, 1> b_tensor = Eigen::TensorMap<const Eigen::Tensor<const T, 1>> (b.data(), b.size());
            const Eigen::Tensor<T, 0> ret = b_tensor.contract(a_tensor, product_dims);
            return ret();
        }

then I get

Code: Select all

/usr/include/eigen3/unsupported/Eigen/CXX11/src/Tensor/TensorContraction.h:418:117: error: no matching function for call to ‘Eigen::internal::general_matrix_vector_product<long int, stan::math::var, Eigen::internal::TensorContractionInputMapper<stan::math::var, long int, 1, Eigen::TensorEvaluator<const Eigen::Tensor<stan::math::var, 1, 0, long int>, Eigen::DefaultDevice>, std::array<long int, 0>, std::array<long int, 1>, 1, true, false, 16>, 0, false, double, Eigen::internal::TensorContractionInputMapper<double, long int, 0, Eigen::TensorEvaluator<const Eigen::Tensor<double, 1>, Eigen::DefaultDevice>, std::array<long int, 0>, std::array<long int, 1>, 2, true, true, 16>, false, 0>::run(const Index&, const Index&, LhsMapper&, RhsMapper&, Eigen::TensorContractionEvaluatorBase<Eigen::TensorEvaluator<const Eigen::TensorContractionOp<const std::array<const Eigen::IndexPair<int>, 1>, const Eigen::Tensor<stan::math::var, 1, 0, long int>, const Eigen::Tensor<double, 1> >, Eigen::DefaultDevice> >::Scalar*&, const Index&, const Scalar&)’
     internal::general_matrix_vector_product<Index,LhsScalar,LhsMapper,ColMajor,false,RhsScalar,RhsMapper,false>::run(
     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^
         rows, cols, lhs, rhs,
         ~~~~~~~~~~~~~~~~~~~~~                                                                                       
         buffer, resIncr, alpha);
         ~~~~~~~~~~~~~~~~~~~~~~~                                                                                     
In file included from /usr/include/eigen3/Eigen/Core:473:0,
                 from /usr/include/eigen3/Eigen/Dense:1,
                 from [redacted]
                 from [redacted]
/usr/include/eigen3/Eigen/src/Core/products/GeneralMatrixVector.h:88:24: note: candidate: static void Eigen::internal::general_matrix_vector_product<Index, LhsScalar, LhsMapper, 0, ConjugateLhs, RhsScalar, RhsMapper, ConjugateRhs, Version>::run(Index, Index, const LhsMapper&, const RhsMapper&, Eigen::internal::general_matrix_vector_product<Index, LhsScalar, LhsMapper, 0, ConjugateLhs, RhsScalar, RhsMapper, ConjugateRhs, Version>::ResScalar*, Index, RhsScalar) [with Index = long int; LhsScalar = stan::math::var; LhsMapper = Eigen::internal::TensorContractionInputMapper<stan::math::var, long int, 1, Eigen::TensorEvaluator<const Eigen::Tensor<stan::math::var, 1, 0, long int>, Eigen::DefaultDevice>, std::array<long int, 0>, std::array<long int, 1>, 1, true, false, 16>; bool ConjugateLhs = false; RhsScalar = double; RhsMapper = Eigen::internal::TensorContractionInputMapper<double, long int, 0, Eigen::TensorEvaluator<const Eigen::Tensor<double, 1>, Eigen::DefaultDevice>, std::array<long int, 0>, std::array<long int, 1>, 2, true, true, 16>; bool ConjugateRhs = false; int Version = 0; Eigen::internal::general_matrix_vector_product<Index, LhsScalar, LhsMapper, 0, ConjugateLhs, RhsScalar, RhsMapper, ConjugateRhs, Version>::ResScalar = stan::math::var]
 EIGEN_DONT_INLINE void general_matrix_vector_product<Index,LhsScalar,LhsMapper,ColMajor,ConjugateLhs,RhsScalar,RhsMapper,ConjugateRhs,Version>::run(
                        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/include/eigen3/Eigen/src/Core/products/GeneralMatrixVector.h:88:24: note:   no known conversion for argument 7 from ‘const Scalar {aka const stan::math::var}’ to ‘double’

Your help would be much appreciated.

Many thanks,

krey