Parallelization of LSE solvers using CUDA

Question

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 )

Answer 1

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.