Tridiagonal Solvers: Parallel Cyclic Reduction
Derivation and implementation of Parallel Cyclic Reduction
Parallel Cyclic Reduction (PCR) (Hockney and Jesshope 1981) is a variation of the Cyclic Reduction algorithm designed for maximum parallelism. While standard Cyclic Reduction (CR) reduces the system size by half at each step and requires a sequential upsweep (back-substitution) phase, PCR eliminates the upsweep phase entirely.
In PCR, very row of the tridiagonal system is updated simultaneously at each step. By the end of \(\log_2 N\) steps, the original tridiagonal matrix is transformed into a diagonal matrix, allowing every \(x_i\) to be solved independently and concurrently. PCR performs significantly more total arithmetic than CR, but it does so in fewer synchronization steps. On some hardware where arithmetic is cheaper than global synchronization, PCR could outperform CR despite the extra math.
This article is intended for readers who are looking for a deep dive into Parallel Cyclic Reduction. For an introduction to the problem, benchmarks, and discussion of algorithm choices, please see the Tridiagonal Solvers article.
3. Linear Algebra Derivation
Consider the tridiagonal equation for any row \(i\):
\[a_i x_{i-1} + b_i x_i + c_i x_{i+1} = d_i \tag{1}\]
At each step \(k\) (where \(distance = 2^k\)), we use rows \(i - distance\) and \(i + distance\) to eliminate the dependencies of row \(i\) on its immediate neighbors.
Let the equations at distance \(h\) be: Eq \(i-h\): \(a_{i-h} x_{i-2h} + b_{i-h} x_{i-h} + c_{i-h} x_i = d_{i-h}\) Eq \(i+h\): \(a_{i+h} x_i + b_{i+h} x_{i+h} + c_{i+h} x_{i+2h} = d_{i+h}\)
By substituting \(x_{i-h}\) and \(x_{i+h}\) into Equation 1 for row \(i\), we derive the new coefficients for row \(i\):
\[\alpha_i = \frac{a_i}{b_{i-h}}, \quad \beta_i = \frac{c_i}{b_{i+h}} \tag{2}\]
\[ \begin{aligned} a'_i &= -\alpha_i a_{i-h} \\ b'_i &= b_i - \alpha_i c_{i-h} - \beta_i a_{i+h} \\ c'_i &= -\beta_i c_{i+h} \\ d'_i &= d_i - \alpha_i d_{i-h} - \beta_i d_{i+h} \end{aligned} \tag{3}\]
After \(\log_2 N\) steps, the \(a'_i\) and \(c'_i\) coefficients become zero (as they refer to indices outside the matrix bounds), leaving: \[b'_i x_i = d'_i \implies x_i = \frac{d'_i}{b'_i} \tag{4}\]
4. Pseudocode for Implementation
function solve_tridiag_pcr(a, b, c, d):
n = length(d)
// Initial setup: pad a and c to length n
curr_a, curr_b, curr_c, curr_d = pad(a, b, c, d)
num_steps = ceil(log2(n))
for k from 0 to num_steps - 1:
distance = 2^k
// This entire block can be executed in parallel for all i
parallel for i from 0 to n - 1:
im = i - distance
ip = i + distance
// Elimination factors
alpha = (im >= 0) ? curr_a[i] / curr_b[im] : 0
beta = (ip < n) ? curr_c[i] / curr_b[ip] : 0
// Update coefficients for current row
new_a[i] = -alpha * (im >= 0 ? curr_a[im] : 0)
new_c[i] = -beta * (ip < n ? curr_c[ip] : 0)
new_b[i] = curr_b[i] - alpha * (im >= 0 ? curr_c[im] : 0)
- beta * (ip < n ? curr_a[ip] : 0)
new_d[i] = curr_d[i] - alpha * (im >= 0 ? curr_d[im] : 0)
- beta * (ip < n ? curr_d[ip] : 0)
curr_a, curr_b, curr_c, curr_d = new_a, new_b, new_c, new_d
return curr_d / curr_b