# Parallel vs Serial Thomas Algorithm

I am currently writing a code that solves a large tridiagonal matrix every iteration and runs for 1,000's of iterations. I am currently using a Thomas algorithm to solve the matrix serially. I found a parallel version of the Thomas algorithm in the book called "Parallel Scientific Computing in C++ and MPI" (which you can find a version pretty easily if you google it).

The thing is when I run the the parallel Thomas taken directly from the book, it actually has a slower run time than the serial algorithm. This seems baffling to me as the whole point of creating a parallel algorithm is to speed it up. I tested both by running the serial and parallel algorithms verbatim from the book and sending the same tridiagonal into the function for each. I am running MPI on a school cluster to run these and submitting a job for the parallel runs to ensure I am using the same processor type for each node. Could it just be the communication time from the layout of the cluster that is causing this? I feel like I must be doing something wrong so any help would be much appreciated. If posting each algorithm would be of use, let me know.

• "Could it just be the communication time from the layout of the cluster that is causing this?" Can you tell us the results of timing your code? Have you measured the running time of serial/communication portions of your parallel code? – Kirill Dec 1 '14 at 23:33
• The run times for just the parallel and serial Thomas algorithm are as follows for a tridiagonal matrix of size 300,000: 10 processors=0.068049, 8 processors=0.057376, 4 processors=0.034412, 2 processors=0.027935, and serially=0.023686. So it appears to increase the run time pretty linearly as the number of processors increases. – mjaisit Dec 2 '14 at 0:07
• I think you should measure computation/communication time separately. Time scaling linearly with number of processors is evidence of it being dominated by communication. Try computing parallel efficiency. – Kirill Dec 2 '14 at 0:14
• Consider locality - each cpu has, for example, an L1 cache - assuming a commodity cpu. Anytime cpu1 cannot find something in its cache it has to ask another cpu to move the data. See: akkadia.org/drepper/cpumemory.pdf – jim mcnamara Dec 3 '14 at 3:38

The Thomas algorithm is very efficient because its operation count is very low and because data accesses are very likely to be cache hits once data is initially read from memory.

There are two loops.

The first loop traverses the data forward. Each element of the lower, main and upper triangle, along with the right-hand-side vector (which is typically also the solution vector upon completion) will be read only once from memory and used several times as part of the computation.

The second loop traverses the data in reverse order, so data that was just used in the first loop is still present in the cache. On most modern microprocessors, you will have a large highest level cache of several MB. This means, even for a system of 300,000 equations, much of the data will still be in cache.

Any parallel tridiagonal solver I know requires more operations than the serial algorithm. If you count the operations in the algorithm you mentioned from Karniadakis and Kirby, you will see it does many more operations, even though they are done in parallel.

In principle, a parallel tridiagonal solver should be faster than the simple serial algorithm, but at problem sizes and processors count much larger than you are using.

You would have a better chance of speeding up your solution if you could solve for your multiple systems simultaneously, assuming some of them are available at the same time. If that is not possible, but your tridiagonal matrix is the same for multiple solves you could save some serial costs by splitting the Thomas algorithm into two parts, factorization and solve.

If none of these situations is possible, you will likely not see speedup from a parallel version of the Thomas algorithm.

• All correct, it's the memory access that kills you. On GPUs I (and others) have published work using "parallel cyclic reduction" (Heller, circa 1980) for solving tridiagonal systems in parallel. It has $O(n log n)$ run time, but is $n$ parallelizable. It works great on GPUs, and on the vector computers of the bygone era, but is quite useless on modern CPUs because of memory issues. – Aurelius Dec 10 '14 at 15:34