# Linear Systems with Multiple Right Hand sides

I am interested in solving a sequence of linear systems of the form: $$A x_i = b_i$$ That is, all the systems use the same matrix $A$ but they have different right hand sides. The matrix $A$ is sparse symmetric positive definite. I am currently using CHOLMOD to solve each system. Obviously, I factorize $A$ once (Cholesky factorization) and then I perform the forward/backward solves sequentially.

I am now trying to parallelize these operations. I would like to find the $x_i$'s in parallel. On my first attempt using pthreads, I couldn't get more than 2X speed-up on a quad-core i7 machine. I am guessing that I couldn't go beyond 2X because the operations that I am trying to parallelize are so fast that their runtime is comparable to the time required to read the Cholseky factorization matrices from the memory and to send them to each core. I am not sure of that as I am not very experienced in parallel programming.

My question is: is there a better way of doing this? Would I expect higher speed-ups if I used other programming paradigms? Also, are you aware of any solvers out there that already perform this parallelization. (I am not interested in packages that parallelize one system solve).

The problem seems embarrassingly parallel and I expected speed-ups very close to 4X on a quad-core machine as there are no dependencies between the tasks.

Thank you

• You didn't say how large your system is but I'm assuming that most of the computational cost is in the factorization of A. How many RHS do you have? I think it would have to be fairly large for you to expect to achieve any reasonable overall speedup by parallelizing this phase of the solution? – Bill Greene Mar 25 '16 at 22:36
• @BillGreene well.. I am not taking the factorization time into account. Assume that this is already done. I am only interested in the speed-up that I can achieve by paralleling the backward/forward solves. So, ideally speaking, the speed-up should be 4X (not counting the factorization time) since there are no dependencies. But of course it will be lower than that in practice, but I didn't expect it to be as low as 2X. The number of right hand sides is in the hundreds by the way and I have up to 10 million variables. – Mohammad Fawaz Mar 25 '16 at 23:00
• There may be a way of accelerating that, but it will require recoding the backward/forward solve routine of CHOLDMOD: if you interleave the bi's and the xi's and write a backward/forward solve routine in AVX, then you will fetch the sparse matrix once only and will benefit from data locality in both the x_i and b_i's (but I'm not 100% sure this will work and it's a lot of effort) – BrunoLevy Mar 26 '16 at 8:39

## 1 Answer

It's possible that your performance is limited by the memory bandwidth of your system. It's not at all uncommon to have a situation where two or three of your cores can use all of the available memory bandwidth and make it impossible to achieve higher parallel speedups because there isn't enough bandwidth to allow all of the cores to work at full speed.

A relatively simple test is to run John McCalpin's STREAM benchmark with 1, 2, 3, and 4 cores to see how performance scales. Don't be surprised if it bogs down after two cores. If your system is bandwidth limited in this way and if your code is memory intensive (typical for sparse direct solvers), then this could explain your results.

You can see some STREAM benchmark results at McCalpin's web site. One result for a generic core-i7-2600 system shows very little improvement in going from one to four cores.

The STREAM benchmark can be found at:

https://www.cs.virginia.edu/stream/