Dear developers,

my name is Riccardo Rossi, i am one of the main developers of the Kratos project (http://www.cimne.com/kratos).

I begin with a disclaimer: I tried first to post on stack overflow, but my post was closed as "too broad" http://stackoverflow.com/questions/3918 ... e-to-eigen

The reason for which i am writing here is that i would like to experiment using eigen within our code basis. So far this is only a personal initiative, and for the moment i would like simply to explore the potential feasibility of switching to eigen.


As of today we are using boost.ublas as our dense linear algebra system and i am currently trying to see if there is any simple way to switch from ublas to eigen without having to touch the whole code basis (which is of the order of 2M lines).

having had a look to the eigen documentation, i realized that the interface is not too different from the one in ublas, so it might be possible to write a wrapper to ease the porting (i appreciate this is not optimal, but porting would be otherwise infeasible for us)

for example if i have (consider MatrixType as a convenient typedef)

MatrixType A,B,C

i currently write in ublas

noalias(A) = prod(B,C);

while i understand that the equivalent of this in eigen would be

A.noalias() = B*C

my idea would be to try to define in the kratos namespace two helper functions

expression_type = noalias( sometype ) ... giving back A.noalias()
expression_type = prod(A,B) giving back the expression A*B

here the big problem is that i do not have a clue of how this function should look like. (particularly which sort of expression it would give back).

this is a little more complicated since in ublas it is valid to write

VectorType x,y;
MatrixType A;
noalias(x) += prod(A,y);

while (as i understand)

x.noalias() += A*y;

would not be valid for ublas

similar functions for which i would need a wrapper are inner_prod outer_prod norm_2 row

a distinct challenge is that ublas matrix has the methods

size1() --> corresponds to eigen's rows()
size2() --> corresponds to cols()
resize(size1, size2, copy_data=true)

my idea here is that to do a quick and dirty test i could hack into eigen's matrix class to add these methods, while at a later stage the only possibility would be to derive my own matrix class from Eigen's matrix and to add the wrapper methods to my derived class. Do you foresee any performance/compatibility issue with this?

other unrelated questions that come to my mind are related to the performance improvement i might get from this. the scenario in which i work is that OPENMP should never be used in the matrix multiplications (it is used one level up, so for me eigen would behave as purely serial). the matrices and vectors i work with are either dynamic of size of the order 20*20 or static of order 3*3, 4*4 (eventually rarely 20*20)

i read in the benchmarks that ublas may work better in this range of sizes, is that so?? is it something expected to change in the near future? (eigen is so much superior to ublas that i might think into changing anyway, but i would like to have a feedback on this)

please tell me very openly of what you think on the porting approach i am considering. also tell me if you believe it is a bad idea overall.

finally i would like to remark that this is a preliminary exploration project, which may well not go anywhere, particularly since i am not too clear if our license (BSD) is compatible with eigen's (MTL) from the point of view of eventual industrial users of our code.

thank you in advance for any feedback you can give Riccardo

For the record, this is a duplicate of: http://stackoverflow.com/questions/3918 ... e-to-eigen