Hello,

I am trying to execute at the fastest possible speed a matrix inversion process (note that in my case I have no other choice than actually fully invert the matrix as I need ALL the X(i,j)).

Note that my Matrix size could go from M=N=5000 to M=N=30000 and is Symmetric/Diagonally dominant. I concluded that it was Positive Semi Definite (PSD).

This is where I am right now:

I am using Eigen's trunk from yesterday.

1) Despite the fact that my matrix is PSD, which makes it factorizable by LDLt method, I opted for the PartialPivLU solver though, as it performs the factorization in parallel (I can see my 4 cores going full blown) while LDLt does not (I did enabled openmp), making it much slower.

2) I call X = solver.solve(eye) where eye = Matrix::Identity(n,n), but, the process is not parallelized. I would expect that, because each column of X is computed independently, the process should be parallelized. But it's not

What about an openmp version of the solve(A) method?


That brings an important general comment:

There are many very interesting and clear summary tables of the different algorithms throughout the Eigen documentation (which is awesome BTW).

But I think, in this same documentation, there is a big lack of information regarding the potential parallelization of the implementation of the algorithm. The fact that the implementation is parallelized or not is quite important IMHO: I've nearly cut by 4 my factorization by switching from LDLt to PartialPivLU!


What do you think?


SR

Hi,

the truth is that currently only general matrix-matrix products are parallelized. Since the cost of PartialPivLU is dominated by such matrix products, this later operation is naturally parallelized, but that's all.

In the near future we plan to parallelize other kind of products, and in particular triangular solves such that LDLT, LLT, and .solve(.) methods benefit from multithreading. However, to get good scaling, this is not as trivial as parallelizing over the columns...

Good idea to add an openmp note in the big solver table.