Cast DMatrix to Eigen Matrix and back

[Peter Listov]

Hello,

Recently I encountered a problem (well documented though) that some of the linear algebra is not so fast in CasADi. That is, for instance, computation of the inverse of 13x13 matrix takes half a second. So I am considering to use Eigen library for time critical sections in a program. To copy data from DM to Eigen I use the Map class:

Example:

DM E = DM::eye(7);
Eigen::Matrix<double, 7, 7> mat = Eigen::Matrix<double, 7, 7>::Map(DM::densify(E).nonzeros().data(), 7, 7);

Would anyone recommend a good way to cast Eigen Matrix back to DM type?

Thanks!

Best regards,

Peter

[Joris Gillis]

Dear Peter,

You should be careful about chaining operations like this. You may be passing a pointer to a temporary object here.
Better save the result of DM::densify(E).nonzeros() first...

Best regards,
    Joris

[Peter Listov]

Indeed. Thanks, Joris!

And the only solution for the reverse problem (Eigen -> DM) I came up with so far is:

1. Copy Eigen Matrix to STL vector of doubles;

2. Call the DM constructor with this vector

3. And then call 'reshape' procedure 

Not sure it's a good one however.

Best regards,

Peter

[Joris Gillis]

An alternative is to allocate a DM::zeros of correct shape, and using the .ptr() member function to copy from Eigen to DM directly.

[Joel Andersson]

Hi!

One side remark: In the case of the "inv" function, note that CasADi's "inv" is using minor expansion, which is very inefficient for larger dimension. The "solve" operation, scales much better:

from casadi import *

A = SX.sym('A', 8, 8);

A_inv1 = inv(A) # very slow

A_inv2 = solve(A, SX.eye(8)) # fast

I created an issue for this since it's a potential pitfall: https://github.com/casadi/casadi/issues/1871

Best regards,

Joel

[Paul Daum]

Hi,

I had the same problem Peter and this seems to be the best way to do this. I also suggest using memcpy for converting DM to Eigen, since it's a little bit faster.
Here's the code I used:

From Eigen to Casadi:

size_t rows = eigen_matrix.rows();
size_t cols = eigen_matrix.cols();

casadi_matrix.resize(rows,cols);
casadi_matrix = casadi::DM::zeros(rows,cols);

std::mempy(casadi_matrix.ptr(), eigen_matrix.data(), sizeof(double)*rows*cols);

From Casadi to Eigen:

size_t rows = casadi_matrix.size1();
size_t cols = casadi_matrix.size2();

eigen_matrix.resize(rows,cols);
eigen_matrix.setZero(rows,cols);

std::memcpy(eigen_matrix.data(), casadi_matrix.ptr(), sizeof(double)*rows*cols);

[Paul Daum]

I just noticed that the way I described for converting a Casadi matrix to an Eigen matrix only works for dense matrices.
If you have a sparse matrix you need to densify it first, as you did in your original post.

[Peter Listov]

Hi,

I've read once that one might lose a bit of Eigen performance when using memcpy due to memory alignment issues. In order to preserve sparsity I use this:

casadi::DM Matrx = ...

casadi::Sparsity SpA = Matrx.get_sparsity();

std::vector<int> output_row, output_col;

SpA.get_triplet(output_row, output_col);

std::vector<double> values = Matrx.get_nonzeros();

using T = Eigen::Triplet<double>;

std::vector<T> TripletList;

TripletList.resize(values.size());

for(int k = 0; k < values.size(); ++k)

    TripletList[k] = T(output_row[k], output_col[k], values[k]);

Eigen::SparseMatrix<double> SpMatrx(Matrx.size1(), Matrx.size2());

SpMatrx.setFromTriplets(TripletList.begin(), TripletList.end());

  

Best regards,

Peter