Tridiagonal Solvers: Cyclic Reduction
Derivation and implementation of Cyclic Reduction
Cyclic Reduction (CR) (Hockney 1965) is a direct method for solving tridiagonal linear systems of equations. It is more suitable for parallel computing architectures as it utilizes simultaneously eliminating odd- or even-indexed unknowns.
Unlike the Thomas algorithm, which is inherently sequential, CR repeatedly transforms the original system into a smaller system of the same form. This reduction process continues until the system becomes small enough to be solved directly, after which the eliminated unknowns are recovered through a back-substitution phase.
This method is most naturally applied when the number of equations \(N\) is of the form \(2^k - 1\) for some integer \(k\). If \(N\) is not of this form, the system might need to be padded with identity rows or handled with special boundary conditions, but the core principle of recursive reduction remains.
This article is intended for readers who are looking for a deep dive into Cyclic Reduction. For an introduction to the problem, benchmarks, and discussion of algorithm choices, please see the Tridiagonal Solvers article.
Derivation of the Algorithm
The reduction process works by eliminating every other unknown (e.g., all even-indexed \(x_i\) if the system is indexed \(0, 1, \dots, N-1\)). This results in a new, smaller tridiagonal system involving only the remaining (odd-indexed) unknowns. This process is then repeated on the smaller system.
For this derivation, consider again our system.
\[ a_i x_{i-1} + b_i x_i + c_i x_{i+1} = d_i \tag{1}\]
with the boundary coefficients \(a_0 = 0\) and \(c_{N-1} = 0\).
where: - \(b_i\) are the elements of the main diagonal. - \(a_i\) are the elements of the sub-diagonal (\(a_0\) is usually defined as 0 or ignored). - \(c_i\) are the elements of the super-diagonal (\(c_{N-1}\) is usually defined as 0 or ignored).
Phase 1: Reduction (Downsweep)
Starting with the general tridiagonal equation Equation 1 and boundary coefficients \(a_0 = 0\) and \(c_{N-1} = 0\), we can manipulate the system to eliminate every other unknown. Consider three consecutive equations where \(i\) is odd:
\[ \begin{aligned} a_{i-1} x_{i-2} + b_{i-1} x_{i-1} + c_{i-1} x_i &= d_{i-1} \\ a_i x_{i-1} + b_i x_i + c_i x_{i+1} &= d_i \\ a_{i+1} x_i + b_{i+1} x_{i+1} + c_{i+1} x_{i+2} &= d_{i+1} \end{aligned} \tag{2}\]
From the first and third equations in Equation 2, we can express the even-indexed unknowns \(x_{i-1}\) and \(x_{i+1}\):
\[ \begin{aligned} x_{i-1} &= (d_{i-1} - a_{i-1} x_{i-2} - c_{i-1} x_i) / b_{i-1} \\ x_{i+1} &= (d_{i+1} - a_{i+1} x_i - c_{i+1} x_{i+2}) / b_{i+1} \end{aligned} \]
Substituting these into the middle equation of Equation 2:
\[ a_i \left( \frac{d_{i-1} - a_{i-1} x_{i-2} - c_{i-1} x_i}{b_{i-1}} \right) + b_i x_i + c_i \left( \frac{d_{i+1} - a_{i+1} x_i - c_{i+1} x_{i+2}}{b_{i+1}} \right) = d_i \]
Rearranging terms to group the remaining odd-indexed unknowns (\(x_{i-2}\), \(x_i\), and \(x_{i+2}\)):
\[ \left( \frac{-a_i a_{i-1}}{b_{i-1}} \right) x_{i-2} + \left( b_i - \frac{a_i c_{i-1}}{b_{i-1}} - \frac{c_i a_{i+1}}{b_{i+1}} \right) x_i + \left( \frac{-c_i c_{i+1}}{b_{i+1}} \right) x_{i+2} = d_i - \frac{a_i d_{i-1}}{b_{i-1}} - \frac{c_i d_{i+1}}{b_{i+1}} \]
This results in a new tridiagonal equation involving only odd-indexed unknowns. We can define the modified coefficients for this reduced system as:
\[ \begin{aligned} a'_i &= \frac{-a_i a_{i-1}}{b_{i-1}} \\ b'_i &= b_i - \frac{a_i c_{i-1}}{b_{i-1}} - \frac{c_i a_{i+1}}{b_{i+1}} \\ c'_i &= \frac{-c_i c_{i+1}}{b_{i+1}} \\ d'_i &= d_i - \frac{a_i d_{i-1}}{b_{i-1}} - \frac{c_i d_{i+1}}{b_{i+1}} \end{aligned} \tag{3}\]
If we define
\[ \begin{aligned} \alpha &= \frac{a_i}{b_{i-1}} \\ \beta &= \frac{c_i}{b_{i+1}} \end{aligned} \tag{4}\]
we can then simplify equation Equation 3 by substituting Equation 4 to get a convenient form:
\[ \begin{aligned} a'_i &= -\alpha a_{i-1} \\ b'_i &= b_i - \alpha c_{i-1} - \beta a_{i+1} \\ c'_i &= -\beta c_{i+1} \\ d'_i &= d_i - \alpha d_{i-1} - \beta d_{i+1} \end{aligned} \tag{5}\]
This reduction step is applied iteratively until a \(1 \times 1\) system remains, which is solved directly.
Phase 2: Substitution (Upsweep)
Once the unknowns of the smallest system are found, we recover the eliminated unknowns in reverse order. For any eliminated even-indexed unknown \(x_j\), we use its original equation:
\[ a_j x_{j-1} + b_j x_j + c_j x_{j+1} = d_j \]
Since \(x_{j-1}\) and \(x_{j+1}\) are either from the reduced system or were recovered in a previous upsweep step, we solve for \(x_j\) directly:
\[ x_j = (d_j - a_j x_{j-1} - c_j x_{j+1}) / b_j \tag{6}\]
Pseudocode for Implementation
To help illustrate the difference between the mathematical logic and the hardware-optimized implementation, we provide both sequential and parallel versions of the Cyclic Reduction algorithm. Both versions assume \(N = 2^k - 1\) for simplicity.
Sequential Pseudocode
The sequential version is useful for understanding the step-by-step transformation of the tridiagonal system.
function solve_tridiag_cyclic_reduction(a, b, c, d):
n = length(d)
// Store coefficients at each level for back-substitution
history = []
curr_a, curr_b, curr_c, curr_d = a, b, c, d // Initialize with original coefficients
// --- Phase 1: Reduction (Downsweep) ---
while length(curr_d) > 1:
history.append((curr_a, curr_b, curr_c, curr_d))
next_a = []
next_b = []
next_c = []
next_d = []
// Iterate over odd indices (1, 3, 5, ...) in the current system to keep them
// for the next level (reduced system)
for i from 1 to length(curr_d) - 2 step 2:
// Calculate factors to eliminate adjacent even-indexed unknowns
alpha_factor = curr_a[i-1] / curr_b[i-1] // Factor for even row i-1
beta_factor = curr_c[i] / curr_b[i+1] // Factor for even row i+1
// New coefficients for the reduced equation (corresponding to original odd-indexed rows)
next_b.append(curr_b[i] - alpha_factor * curr_c[i-1] - beta_factor * curr_a[i])
next_d.append(curr_d[i] - alpha_factor * curr_d[i-1] - beta_factor * curr_d[i+1])
// New sub-diagonal element for the next level
if i > 1:
next_a.append(-alpha_factor * curr_a[i-2])
// New super-diagonal element for the next level
if i < length(curr_d) - 2:
next_c.append(-beta_factor * curr_c[i+1])
curr_a, curr_b, curr_c, curr_d = next_a, next_b, next_c, next_d
// Solve the final 1x1 system
x = [curr_d[0] / curr_b[0]]
// --- Phase 2: Substitution (Upsweep) ---
for prev_a, prev_b, prev_c, prev_d in reversed(history):
new_x = array_of_zeros(length(prev_d))
// Fill in the known values from the level above (which are odd-indexed here)
for j from 0 to length(x) - 1:
new_x[2*j + 1] = x[j]
// Solve for the "missing" even-indexed values
for i from 0 to length(prev_d) - 1 step 2: // Iterate over even indices
rhs = prev_d[i]
if i > 0: // Check for sub-diagonal term (x_{i-1})
rhs -= prev_a[i-1] * new_x[i-1]
if i < length(prev_d) - 1: // Check for super-diagonal term (x_{i+1})
rhs -= prev_c[i] * new_x[i+1]
new_x[i] = rhs / prev_b[i]
x = new_x // Update x for the next upsweep level
return xParallel Pseudocode
The parallel version highlights the algorithm’s inherent efficiency on modern hardware by using parallel for constructs for both the reduction and recovery phases.
function solve_tridiag_parallel_cyclic_reduction(a, b, c, d):
n = length(d)
history = []
curr_a, curr_b, curr_c, curr_d = a, b, c, d
// --- Phase 1: Reduction (Downsweep) ---
while length(curr_d) > 1:
history.append((curr_a, curr_b, curr_c, curr_d))
// Pre-allocate space for the reduced system
m = (length(curr_d) - 1) / 2
next_a = array_of_zeros(m - 1)
next_b = array_of_zeros(m)
next_c = array_of_zeros(m - 1)
next_d = array_of_zeros(m)
// All reduction steps can be performed in parallel
parallel for i from 1 to length(curr_d) - 2 step 2:
j = (i - 1) / 2 // Index in the reduced system
alpha = curr_a[i-1] / curr_b[i-1]
beta = curr_c[i] / curr_b[i+1]
// Compute new coefficients concurrently
next_b[j] = curr_b[i] - alpha * curr_c[i-1] - beta * curr_a[i]
next_d[j] = curr_d[i] - alpha * curr_d[i-1] - beta * curr_d[i+1]
if i > 1:
next_a[j-1] = -alpha * curr_a[i-2]
if i < length(curr_d) - 2:
next_c[j] = -beta * curr_c[i+1]
curr_a, curr_b, curr_c, curr_d = next_a, next_b, next_c, next_d
// Solve the final 1x1 system
x = [curr_d[0] / curr_b[0]]
// --- Phase 2: Substitution (Upsweep) ---
for prev_a, prev_b, prev_c, prev_d in reversed(history):
new_x = array_of_zeros(length(prev_d))
// Populate known odd-indexed values in parallel
parallel for j from 0 to length(x) - 1:
new_x[2*j + 1] = x[j]
// Concurrent recovery of all even-indexed values
parallel for i from 0 to length(prev_d) - 1 step 2:
rhs = prev_d[i]
if i > 0:
rhs -= prev_a[i-1] * new_x[i-1]
if i < length(prev_d) - 1:
rhs -= prev_c[i] * new_x[i+1]
new_x[i] = rhs / prev_b[i]
x = new_x
return x