We are currently evaluating Eigen as a replacement for an in-house Linear Algebra library.

We have many functions that take an "abstract" array of a given data type and dimension (1D/2D) as parameters. Users can then pass either concrete arrays holding their own data, or views on those arrays, similar to Eigen blocks. For example, a function taking a vector can also take the row or column of a matrix.

Since the code base is organized in class hierarchies using virtual functions, we cannot implement those functions as templates. Is there a way to express "any matrix with direct data access" as a data type in Eigen?

Example pseudo-code:

Code: Select all

// Is there a type in Eigen I can substitute for SomeDynamicMatrix that
// allows both Eigen::Matrix<...> and Eigen::Block<...>, without causing
// copies?
virtual int func(const SomeDynamicMatrix<float>& input, SomeDynamicMatrix<float>& output);


see this page: http://eigen.tuxfamily.org/dox/TopicFun ... Types.html

Thank you for your reply.

I have read the page you are referring me to, but my question still remains: If I want to use ordinary functions, not template functions, I cannot use MatrixBase. If I take a Matrix&, users cannot pass blocks. If I take a Block, users cannot pass matrices (without taking a block spanning the whole matrix) and I suspect that "special" blocks won't work either. Does that mean that there is no built-in way?

I'm thinking of a kind of "array descriptor", holding a data pointer and stride information, that can be constructed from arrays, matrices and blocks.

One solution I can think of is to create a type (let's call it View) that takes an Eigen expression and holds a Map<MatrixXt>, created using data(), rows() and cols(). This could than be used as follows:

Code: Select all

void func(const ViewNt& input, ViewNt& output)
{
    // Access Eigen API by converting the views to Eigen::Map
    output.map() = input.map();
}


Is there a problem with that approach?

This is indeed something we would like to add soon. Basically it boils dow to adding appropriate ctors to Map<>.

I talked about that in this topic : viewtopic.php?f=74&t=96780

In my current project I exactly use the approach you described (a View class deriving from Map + appropriate ctor), for passing matrices blocks to non-template functions & storing part of vectors into containers without copying the original content.

You'd be aware that with this approach you can sometimes loose performances.

I wrote this small benchmark to illustrate my words :

Code: Select all

MatrixXd mat1(1000,1000);
MatrixXd mat2(1000,1000);
MatrixXd mat3(1000,1000);

// To make sure there is no NaN !
mat1.setRandom();
mat2.setRandom();
mat3.setRandom();

// View is basically a Map with dynamic strides + special ctor
View<VectorXd> view1(mat1.col(50));
View<VectorXd> view2(mat2.col(50));
View<VectorXd> view3(mat3.col(50));

clock_t start = clock();

for(int i = 0 ; i < 100000 ;i++) {

   // Direct multiplication (A)
   mat3.col(50) = mat1.col(50)*mat2.col(50);

}

printf("A = %dms\n", 1000*(clock()-start)/CLOCKS_PER_SEC);

start = clock();

for(int i = 0 ; i < 100000 ;i++) {

   // Multiplication using views (B)
   view3 = view1*view2;
}

printf("B = %dms\n", 1000*(clock()-start)/CLOCKS_PER_SEC);


Compiler : MSVC2010 /Release.
CPU : Intel Core i7-720QM 1.6Ghz

Output :
A = 839ms
B = 1578ms

> As you can see, there is a huge difference(~2x) between direct multiplication & multiplication of maps. Maybe it's because the strides are not known at compile time, but it's just my hypothesis.

And if I replace col() by row() :
A = 2803ms
B = 3006ms

> Relatively small difference, because the View's cost is certainly hidden by the cache misses (the matrices are col major ordered by default, so a row is not a continuous block a memory)

Conclusion :

I made this error when I started the dev of my application : I created a View class and then I was very happy because template functions were not needed anymore. I know that your case is different because of virtuals, but you should think about this small benchmark if performances are critical.

Maybe you could use a Matrix& when you're sure that the user will pass a plain matrix for example.

This is indeed something we would like to add soon. Basically it boils dow to adding appropriate ctors to Map<>.

Since you seem confident that adding this functionality to Map<> is rather simple, can we expect that this will be added in the near future?

if you declared your View<> to be as generic as possible, i.e., with a dynamic inner stride, then Eigen cannot vectorize anymore, whence the slowdown. Perhaps enforcing an inner stride of 1, would be a good compromise.

fbeyer wrote:Since you seem confident that adding this functionality to Map<> is rather simple, can we expect that this will be added in the near future?

that depends how you define "near"... we have very little time these days.... But if someone come up with a clean patch with thorough unit tests....

ggael wrote:But if someone come up with a clean patch with thorough unit tests....

I understand. I see what I can do but I'm afraid I've too little knowledge of Eigen's "black magic" (yet?).

For the record, the Ref<> class covers this use case: http://eigen.tuxfamily.org/dox/classEigen_1_1Ref.html