Compilation error when multiplying matrix of ceres::Jet by matrix of double
Summary
Hi,
I am using Eigen for matrix multiplication in Ceres solver.
Suppose I have two matrices, one consists of ceres::Jet<double, N>
, whereas the other consists of double
. When I try to calculate the product of the two matrices, I'm getting a compilation error. In short, the error boils down to an attempt of Jet-to-double conversion.
This question has originally been posted on Ceres Solver Group: https://groups.google.com/g/ceres-solver/c/bMH73wZ0hoE/m/EChU22cGAgAJ
It's rather a Ceres-Eigen compatibility issue, not a bug report. Correct me, if this issue should be changed to a feature request.
Environment
- Operating System : Linux
- Architecture : x64
- Eigen Version : 3.4.0
- Ceres Solver Version : 2.0.0
- Compiler Version : LLVM 13.0.0
Minimal Example
Here's an example of code that fails to compile:
typedef ceres::Jet<double, 2> J;
Eigen::MatrixX<J> A;
A.resize(2, 2);
A << (J)1., (J)1., (J)1., (J)1.;
Eigen::MatrixXd B;
B.resize(2, 2);
B << 1., 1., 1., 1.;
Eigen::MatrixX<J> C;
C.noalias() = A * B;
Note, that in this simple case, I could use Matrix2 and Matrix2d instead, and it would work fine. However, I need to use matrices of much higher size, so it's not possible to use compilation-time defined size.
What's more, if both of matrices consist of Jets, the example compiles fine. However, for performance reason, I don't want to convert the matrix of doubles to Jets just to make it compile.
It's also worth noting, that the multiplication works, if matrix of Jets is a vector. For now, the workaround I'm using is to just loop over the rows like that:
C.resize(A.rows(), B.cols());
for (int i = 0; i < A.rows(); i++) C.row(i).noalias() = A.row(i) * B;
It's far from perfect though.
Ceres Solver community has already been working on integration with Eigen. However, it's not possible to fix the described issue entirely on the side of Ceres Solver. As Sergiu Deitsch mentioned (see https://groups.google.com/g/ceres-solver/c/bMH73wZ0hoE ) - "It looks like dynamic size matrices use a variant of GEMM which requires partially specializing Eigen::internal::gebp_traits and the corresponding packet math for mixed type operations. [...] Given the necessary interfaces live in an internal namespace [...], I don't believe it is currently feasible to provide a workaround in Ceres."
Does anyone know how to solve that issue? Thanks in advance
Error Messages
A fragment of compilation output of the first example looks like this:
include/Eigen/src/Core/util/BlasUtil.h:43:79: error: cannot convert 'const ceres::Jet<double, 2>' to 'double' without a conversion operator
EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE To run(const From& x) { return To(x); }
include/Eigen/src/Core/GeneralProduct.h:246:66: note: in instantiation of member function 'Eigen::internal::get_factor<ceres::Jet<double, 2>, double>::run' requested here
RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha);
[...]
include/Eigen/src/Core/products/GeneralMatrixMatrix.h:445:7: note: in instantiation of function template specialization 'Eigen::internal::generic_product_impl<Eigen::Matrix<ceres::Jet<double, 2>, -1, -1, 0>, Eigen::Matrix<double, -1, -1, 0>, Eigen::DenseShape, Eigen::DenseShape, 8>::scaleAndAddTo<Eigen::Matrix<ceres::Jet<double, 2>, -1, -1, 0>>' requested here
scaleAndAddTo(dst, lhs, rhs, Scalar(1));
[...]
main.cpp:24:15: note: in instantiation of function template specialization 'Eigen::NoAlias<Eigen::Matrix<ceres::Jet<double, 2>, -1, -1, 0>, Eigen::MatrixBase>::operator=<Eigen::Product<Eigen::Matrix<ceres::Jet<double, 2>, -1, -1, 0>, Eigen::Matrix<double, -1, -1, 0>, 0>>' requested here
[build] C.noalias() = A * B;