Re: [eigen] triangular solve

[ Thread Index | Date Index | More Archives ]

On 21.12.2012 23:31, Gael Guennebaud wrote:
Since your matrices are quite large, the cost of this temporary is
completely negligible, you can still avoid one by doing:

MatrixXf tmp(this->matrixL());
tmp = total_cov.matrixL().solve(tmp); // in-place solving
float denom = tmp.squaredNorm();

Higher performance could be achieved by taking into account that tmp is
triangular and so total_cov.matrixL().solve(tmp) is triangular too.
However, Eigen is not able to do so yet, and I'm not aware of any library
being able to do so.

You can try something like this (pseudocode/not tested):

MatrixXf tmp(this->matrixL());
for([i0, i1]; /* seperate [0..n] into intervals */) {
total_cov.matrixL().bottomRightCorner(n-i1, n-i1).solveInPlace(tmp.block(i0, i0, n-i0, i1-i0));

maybe somewhere a .triangularView() has to be inserted. If you go through the matrix column by column, you can even save the copy. The question would be how that impacts vectorization -- I guess we need better/"smarter" support for triangular matrices, eventually. E.g. saving big triangular matrices wastes a lot of space at the moment. Another nice feature would be things such as multiplyInPlace(), inverse()/invertInPlace() (for small matrices it can be worthwile to invert them once and multiply, thus saving high-cost divisions).


Dipl.-Inf. Christoph Hertzberg
Cartesium 0.049
Universität Bremen
Enrique-Schmidt-Straße 5
28359 Bremen

Tel: +49 (421) 218-64252

Mail converted by MHonArc 2.6.19+