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

Question

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.

Logogriph

Answer 1

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();

Lotz

Your first solution doesn't compile (ISO C++ forbids declaration of ‘type name’ with no type [-fpermissive). Your second solution works, but could you clarify please if the copy made is shallow or deep? I also need to apply the same procedure to matrix c after it has been unconst-ed (see edit, please) – Logogriph

In both solution, if the storage of the input expression matches your needs, then actual_b will just be a reference, otherwise the expression will be evaluated (deep copy). – Lotz

Thanks, it would seem I can't do better than that. – Logogriph

I am getting this error when I try this solution: no known conversion for argument 1 from ‘Eigen::MatrixBase<Eigen::Map<Eigen::Matrix<double, -1, -1, 1, -1, -1>, 0, Eigen::Stride<0, 0> > >’ to ‘Eigen::Ref<Eigen::Matrix<double, -1, -1> >’. Seems Eigen can't convert to a Ref all the time? – Pages

You are trying to reference a row-major mapped matrix by a column-major ref, this cannot work without a deep copy and thus a Ref<const MatrixXd>. – Lotz

Answer 2

I guess the issue is this: you are not allowed to get a pointer to data of a MatrixBase<Derived>, since the latter can be any kind of expression in Eigen, like a product of matrices for example. To get a pointer you probably have to first implicitly convert the MatrixBase<Derived> into a Matrix<Scalar, Dynamic, Dynamic>, then use the data() member of the latter.

So you can create a deep copy of the expression, i.e. use something like

Eigen::Matrix<typename Derived::Scalar, Eigen::Dynamic, Eigen::Dynamic tmp = b;

then use

tmp.data()

This code works now

#include <Eigen/Dense>
#include <iostream>

template<typename Derived>
void use_data\
(const Eigen::MatrixBase<Derived>& mat)
{

    Eigen::Matrix<typename Derived::Scalar, Eigen::Dynamic, Eigen::Dynamic>tmp = mat();
    typename Derived::Scalar* p = tmp.data();

    std::cout << std::endl;
    for(std::size_t i = 0; i < tmp.size(); i++)
        std::cout << *(p+i) << " ";

}

int main()
{
    Eigen::MatrixXd A = Eigen::MatrixXd::Random(2, 2); 
    Eigen::MatrixXd B = Eigen::MatrixXd::Random(2, 2); 

    // now A*B is an expression, of type MatrixBase<EigenSum....>
    use_data(A + B);
}

Bucentaur

Your solution doesn't seem to work for me... float * bMat = (static_cast<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic> >(b)).data(); It segfaults the first time I try to access any elements of the array. – Logogriph

Yeah I guess it's not the way to go... What happens if you cout the matrix before getting the data pointer? – Bucentaur

Segfaults... I guess I should've pasted it: float * bMat = (static_cast<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic> >(b)).data(); std::cout << "Bmat 0 is: "<< std::endl; std::cout << bMat[0] << std::endl; Segfaults as soon as I try to access it for the cout – Logogriph

@XapaJIaMnu, no, I was asking what if you cout << b; before getting the data pointer from it. – Bucentaur

makes no difference, still segfaults. – Logogriph

That's pretty much what I am doing at the moment (check the second code sample), but this makes an additional memory copy, which is very slow, especially when concerning very large matrices. I am looking for a solution that would allow me to access the array without that. – Logogriph

Answer 3

There are an easy solution to solve your question, combine EigenMap, &a(0, 0) and const_cast you could resue the buffer of the MatrixBase.

Example :

template<typename Derived1,
         typename Derived2>
void example(Eigen::MatrixBase<Derived1> const &input,
             Eigen::MatrixBase<Derived2> const &output)
{
    static_assert(std::is_same<Derived1::Scalar, Derived2::Scalar>::value,
                  "Data type of matrix input, weight, bias and output should be the same");   

        using Scalar = typename Derived3::Scalar;
        using MatType = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
        using Mapper = Eigen::Map<const MatType, Eigen::Aligned>;

        //in the worst case, you can do const_cast<Scalar *> on
        //&bias(0, 0).That is, if you cannot explicitly define the Map
        //type as const        
        Mapper Map(&input(0, 0), input.size());
        output.colwise() += Map; 
    }
}

I try it on windows 8, vc2013 32bits, Eigen version is 3.2.5, no segmentation fault occur(yet), every things looks perfectly fine. I also check the address of the Map, it is same as the original input. You can verify it with another example

#include <Eigen/Dense>

#include <iostream>

template<typename Derived>
void example_2(Eigen::MatrixBase<Derived> &input)
{
    using Scalar = decltype(input[0]);

    Eigen::Map<Derived> map(&input(0, 0),
                            input.rows(),
                            input.cols());
    map(0, 0) = 300;
}

int main()
{           
    Eigen::MatrixXd mat(2, 2);
    mat<<0, 1, 2, 3;
    example_2(mat);
    std::cout<<mat<<"\n\n";    

    return 0;
}

The first element of mat will be "300"

Lontson

The problem in my case was that MatrixBase<Derived> was lazily evaluated, in which case the underlying array didn't exist at the time I was trying to access it. – Logogriph