Eigen, how to access the underlying array of a MatrixBase<Derived>

XapaJIaMnu picture XapaJIaMnu · Aug 2, 2014 · Viewed 8.5k times · Source

I need to access the array that contains the data of a MatrixBase Eigen matrix.

The Eigen library has the data() method which returns a pointer to an array, however it is only accessible from a Matrix type. The MatrixBase doesn't have a similar method, even though the MatrixBase class is supposed to act as a template and the actual type should be just a Matrix. If I try to access MatrixBase.data() I get a compile time error:

template <typename ScalarA, typename Index, typename DerivedB, typename DerivedC>
void uscgemv(float alpha, 
     const USCMatrix<ScalarA,Index> &a,
     const MatrixBase<DerivedB> &b,
     const MatrixBase<DerivedC> &c_const)
{
    //...some code
    float * bMat = b.data();
    ///more code
}

This code produces the following compile time error.

error: ‘const class Eigen::MatrixBase<Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<float>, Eigen::Matrix<float, -1, 1> > >’ has no member named ‘data’
 float * bMat = b.data();

So I have to resort to gimmicks such as...

float * bMat;
int bRows = b.rows();
int bCols = b.cols();
mallocPinnedMemory(&bMat, bRows*bCols*sizeof(float));
Eigen::Map<Matrix<float, Dynamic, Dynamic> > bmat_temp(bMat, bRows, bCols);
bmat_temp = b;  //THis is SLOW, we should avoid it.

Then I can access the bMat array...

Those copies back-and-forth are the biggest cost in the gpu matrix multiplication, as I essentially I have to make an extra copy, before even coping to the device...

I can't use Eigen-magma, as this is sparse matrix-in-a-weird-format to a dense matrix (or sometimes vector) multiplication so I can't use any of the automatic gpu functions there. Also I would much rather not declare the matrices as something else, because that would require changing A LOT of lines of code across the whole program (which I didn't write).

EDIT: A static cast solution was proposed:

float * bMat = (static_cast<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic> >(b)).data();

However I get segfault the first time I try to access an element of the array bMat.

EDIT 2: I'm looking for a zero copy way to access the underlying arrays. I need to only be able to read b, but I also need to able to write to c. Currently c is unconst-d with the following macro:

#define UNCONST(t,c,uc) Eigen::MatrixBase<t> &uc = const_cast<Eigen::MatrixBase<t>&>(c);

EDIT 3: After cross posting to Eigen Forums it would seem I can't do better than the suggested answer.

Answer

ggael picture ggael · Aug 2, 2014

MatrixBase is the base class of any dense expression. It does not necessarily correspond to an object with storage. For instance, can be the abstract representation of A+B, or in your case the abstract representation of a vector with constant values. You can make uscgemv accepts only expression having appropriate storage using the Ref<> class, e.g.:

template <typename ScalarA, typename Index>
void uscgemv(float alpha, 
 const USCMatrix<ScalarA,Index> &a,
 Ref<const VectorXf> b,
 Ref<VectorXf> c);

If the third argument does not match the storage of a VectorXf then it will be evaluated for you. Then you can safely call b.data(). To keep the scalar type of b generic, you can still declare it as MatrixBase<DerivedB>& and then copy it into a Ref<const Matrix<typename DerivedB::Scalar, DerivedB::RowsAtCompileTime, DerivedB::ColsAtCompileTime> >:

typedef Ref<const Matrix<typename DerivedB::Scalar,  DerivedB::RowsAtCompileTime, DerivedB::ColsAtCompileTime> > RefB;
RefB actual_b(b);
actual_b.data();