Re: [eigen] triangular solve

[Christoph Hertzberg <chtz@xxxxxxxxxxxxxxxxxxxxxxxx>]

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).


Christoph


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

Tel: +49 (421) 218-64252
----------------------------------------------