Understanding solveInPlace operation in Eigen

I was trying to explore the option of "solveInPlace()" function while using LLT in Eigen3.3.7 to speed up the matrix inverse computation in my application. I used the following code to test it.

    int main()
    {
        const int M=3;

        Eigen::Matrix<MyType,Eigen::Dynamic,Eigen::Dynamic> R = Eigen::Matrix<MyType,Eigen::Dynamic,Eigen::Dynamic>::Zero(M,M);
        // to make sure full rank
        for(int i=0; i<M*2; i++)
        {
            const Eigen::Matrix<MyType, Eigen::Dynamic,1> tmp = Eigen::Matrix<MyType,Eigen::Dynamic,1>::Random(M);
            R += tmp*tmp.transpose();
        }

        std::cout<<"R \n";
        std::cout<<R<<std::endl;
        decltype (R) R0 =  R; // saving for later comparison



        Eigen::LLT<Eigen::Ref<Eigen::Matrix<MyType,Eigen::Dynamic,Eigen::Dynamic> > > myllt(R);
        const Eigen::Matrix<MyType,Eigen::Dynamic,Eigen::Dynamic> I = Eigen::Matrix<MyType,Eigen::Dynamic,Eigen::Dynamic>::Identity(R.rows(), R.cols());

        myllt.solveInPlace(I);

        std::cout<<"I: "<<I<<std::endl;
        std::cout<<"Prod InPlace: \n"<<R0*I<<std::endl;


        return 0;
}

After reading the Eigen documentation, I thought that the input matrix (here "R") will be modified while computing the transform. To my surprise, I found that the results is store in "I". This was not expected as I defined "I" as a constant. Please provide an explanation for this behaviour.

Solution

The simple non-compiler answer would be that you're asking for the LLT to solve in-place (i.e. in the passed parameter) so what would you expect the result to be? Apparently, you would expect it to be a compiler error, as the "in-place" means change the parameter, but you're passing a const object.

So, if we search the Eigen docs for solveInPlace, we find the only item that takes a const reference to have the following note:

"in-place" version of TriangularView::solve() where the result is written in other

Warning
The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here. This function will const_cast it, so constness isn't honored here.

The non-in-place option would be:

R = myllt.solve(I);

but that won't really speed up the calculation. In any case, benchmark before you decide that you need the in-place option.

You're question is in place, as what const_cast is meant to do is strip references/pointers of their const-ness iff the underlying variable is not const qualified* (cppref). If you were to write some examples

const int i = 4;
int& iRef = const_cast<int&>(i); // UB, i is actually const
std::cout << i; // Prints "I want coffee", or it can as we like UB
int j = 4;
const int& jRef = j;
const_cast<int&>(jRef)++; // Legal. Underlying variable is not const.
std::cout << j; // Prints 5

The case with i may well work as expected or not, we're dependent on each implementation/compiler. It may work with gcc but not with clang or MSVC. There are no guarantees. As you are indirectly invoking UB in your example, the compiler can choose to do what you expect or something else entirely.

*Technically it's the modification that's UB, not the const_cast itself.