# Parallelization of LSE solvers using CUDA

I want to know methods which are fully parallelizable on CUDA architecture. I have implemented the Jacobi and Conjugate Gradient methods and now Im thinking about the Bi-Conjugate gradient method. I was using the cublas library to simplify some basic algebra operations. Someone said me that Gauss-Seidel is semi-parallelizable. I am new on CUDA, so I cannot understand some issues in the Cholesky descomposition or Gaussian-Elimination papers like:

Gaussian Elimination

Cholesky Descomposition

I can't understand the new algorithms, so I am saying that I have to use the classic numerical algorithms(I am not able to "create" that kind of algorithms), something like.

The point is, does there exist any "theory" about parallelizing of those algorithms? Because it was easy for me the Jacobi parallelization and conjugate gradient methods, there was some saxpy and gemv operations inside the algorithms.

For example in the Cholesky's case the inner loop depends of the first, how we can solve this?

for k = 1 to n
...
for j = k + 1 to n

PS: I always developed the algorithms on C++ and from there I started to parallelizing process( it means the process of "discovering" how to parallelize )

• Hello to scicomp! I have changed the title of your question and I hope others agree that it suits better to your question. – shuhalo Sep 19 '12 at 10:02
• I posted just a simple, common sense answer to your question, which is very broad in scope. If you are interested in more details you should restrict your area of interest. – Stefano M Sep 19 '12 at 11:16

If I understand your question correctly you are interested in the GPU implementation of direct methods for the solution of dense linear systems of equations.

Let me first say that implementing efficient, accurate, and robust linear algebra routines requires considerable expertise, so your first option should be using an existing library, e.g. MAGMA.

If you are interested in your own implementation for didactical purposes, I would suggest you to focus on a model problem, e.g. Cholesky factorization ($A = C^T C$ where $A$ is symmetric positive definite and $C$ is upper triangular), and study the data dependencies between the $c_{ij}$ elements of $C$.

You will soon learn that a "full" parallelization is impossible, but that are several ways to lay out the computation exploiting the GPU parallelism. Googling for "GPU parallel dense cholesky" will show you a vast literature on the argument, so you should not pretend to catch up in a small amount of time.

A few suggestions:

• remember that there are several algorithmic variants of the Cholesky factorization, a good book on the subject should show you more than one,
• do not focus only on the computations but also on the memory layout to reduce latencies
• pay attention on how you implementation scales for increasing matrix sizes
• start by comparing simple variants, and do not try to optimize from the very first try
• don't be disappointed if your implementation is much less efficient than an established library, but rather study their code to learn.

As what regards CUDA, it adopts a programming model that is tailored to a very specific architecture, with less abstraction than other frameworks like OpenCL or API's like MPI. Programming has to be "hardware aware" to obtain good efficiency, since there is no intermediate layer, doing optimization for you. Again memory and arithmetic latency hiding is the key factor, so there is no simple theoretical result or "parallelization model" to invoke, as far as I know. The "CUDA C Best Practices Guide" is a good (and almost necessary) starting point.

• Thank you. I'm developing a LSE library, the problem is I'm starting to develop with CUDA. A question: parallelization on CUDA is a little different between others technologies right? Is "parallelization theory" on CUDA? or Should have to learn CUDA specifics? – facunvd Sep 19 '12 at 11:57
• @facunvd see my revised post. – Stefano M Sep 19 '12 at 12:50
• I forgot to explain why I'm suggesting Cholesky: it is numerically stable without pivoting, so you will not have to deal with the difficulties linked to $LU$ factorization. – Stefano M Sep 19 '12 at 13:09
• Methods which use pivoting need "researching"? – facunvd Sep 19 '12 at 14:25
• @facunvd pivoting is not so straightforward to implement as Cholesky. – Stefano M Sep 19 '12 at 14:35