# 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