Efficient methods for solving thousands of right-hand-sides?
 

Hello, in my application I have to solve a large sparse SPD linear system for thousands of right-hand-sides as part of a numerical simulation package. Currently I am using matlab, and the algorithm uses metis pre-ordering, cholesky factorization, and then uses parfor loops to parallelize the forward/back substitution solves to as many threads as there are CPU cores. The implementation is surprisingly efficient, but there are several drawbacks, most notably that the matlab parfor creates a copy of the sparse factorization of the system to each thread, which is a large overhead as the matrix factorization is many gigabytes.

I am coding in C++ now, and exploring tools and libraries for simultaneously solving thousands of right-hand-sides on a standard multi-core desktop PC (i7 quadcore type machine).

My current approach is going to be to use metis+cholmod to factor the system, and then use OpenMP to possibly parallelize the solution, but I would welcome any other suggestions. For now I am trying to keep the implementation to be non-mpi, that is, running on a single machine, but will likely explore using MPI in the future.

Cheers,

Bob

I don't know what SPD means, can't you just use LU decomposition?
That's very efficient for many RHsides

How are the thousands of solutions going to be used ? Perhaps there is a better way to do it if you explain your situation.

i have written smoothed multigrid lib in c++ that i use in cfd code. I intend to put it online as open source.
(This is based on algo as explained in pyAMG code which is already available as opensource).