c++ - Best way to solve large block diagonal system with Eigen?

Question

I have got this problem: a large 10000x10000 block diagonal matrix M with square blocks of various sizes (for example 1000 blocks, all of size 10), and I want to solve a linear system using the Eigen library.

I think the best way to do this is to loop over each M.block(i, j, p, q) store it into a dense matrix, create a vector of Eigen solvers and compute the decomposition. Then solve the systems in parallel with openmp maybe.

How does one go about doing this, the only way is copying all the .block(i, j, p, q) into dense matrices? then using a vector of Eigen solvers? Can I just store the decomposition into the solvers?

Is there an automatic way to do this in Eigen? Will solving the large sparse linear system do the same (in the sense that it will analyze the pattern and use a built in block diagonal solver?).

Thanks for clarifying!

PS: also I know eigen uses parallelization if the program is compiled with -fopenmp, but this would mean that each 10x10 system is solver in paraller right? Can I switch parallelization off and do it myself on a global system wide basis?

question from:https://stackoverflow.com/questions/65932524/best-way-to-solve-large-block-diagonal-system-with-eigen

Answer 1

Waitting for answers