2D Parallel Approach in Structured Grid Algorithms

• The paper presents a novel 2D parallel approach that preserves convergence by using overlapping domain decomposition and alternating sweep directions.
• Two-dimensional parallel approaches are computational frameworks that leverage dual-axis parallelism in structured grids, enhancing data locality and cache efficiency.
• The method employs localized interface synchronization and matrix splitting to achieve iteration counts similar to sequential sweeps while enabling scalable multicore performance.

A two-dimensional parallel approach refers to algorithms and computational frameworks that exploit parallelism along two distinct, often complementary, axes within a computational or physical system. In contemporary scientific computing and engineering, this concept frequently arises in the context of structured grids, domain decomposition, multi-dimensional physical simulation, and high-performance code optimization—particularly where one-dimensional parallel techniques reach scalability or efficiency limits. In the mathematical and algorithmic literature, two-dimensional parallel methods are rigorously constructed to maintain numerical stability and convergence rates, while offering significant speedups and improved memory/cache behavior relative to purely sequential or naively parallelized alternatives.

1. Foundations: Sequential Sweeping and the Need for Parallelization

Numerous numerical algorithms on structured grids—such as preconditioning via SOR (Successive Over-Relaxation), ILU (Incomplete LU) smoothing, and solution of the eikonal equation using fast sweeping—rely inherently on ordered, sequential updates ("sequential sweeping") dictated by grid topology. These methods, optimized for single-core execution, yield rapid convergence when computational dependencies are locally satisfied. However, their strictly sequential update order fundamentally resists parallelization: naïve attempts to distribute the workload typically degrade convergence, introduce decoupling of critical dependencies, and reduce overall efficiency (Tavakoli, 2010).

The impetus for developing two-dimensional parallel approaches arises from the demands of modern large-scale simulation, where exploiting multicore and distributed architectures is essential. In two-dimensional domains (e.g., rectangles discretized into regular grids), the challenge is to orchestrate concurrent computation such that local sweep dependencies, interface couplings, and global convergence are preserved—necessitating nontrivial algorithmic innovation.

2. Methodological Core: Overlapping Domain Decomposition and Multi-Frontal Sweeping

The paradigm outlined in (Tavakoli, 2010) decomposes the global two-dimensional domain Ω ⊂ ℝ² into a Cartesian array of sub-domains (e.g., a pₓ×p_y tiling). Within each sub-domain, a local sweeping procedure is performed in a prescribed direction. Crucially, the direction of sweeping is alternated between iterations—e.g., alternating between NE→SW and NW→SE diagonals—ensuring that, at each step, newly computed boundary/interface values are preferentially incorporated from the most recent updates of neighboring sub-domains.

To uphold consistency and tightly couple adjacent regions, the method enforces local synchronization at sub-domain interfaces. Specifically, small systems of equations (e.g., 2×2 for edges, 4×4 for corners shared by sub-domains) are solved to update interface values using available data from both internal and neighboring sub-domain sweeps. After this interface update, traditional SOR or Gauss-Seidel sweeps proceed on the (now partially decoupled) interiors.

The algorithm can be formalized as:

1. Partition Ω → {Ω₁, ..., Ω_pₓp_y}
2. For each sub-domain Ω_k:
   • Sweep in assigned direction, updating internal nodes
   • At interface nodes, locally solve small linear systems using neighboring (potentially newly updated) values
3. Alternate sweeping directions for successive global iterations

This hybrid of overlapping domain decomposition and multi-frontal partial coupling constitutes the main methodological advance, preserving strong local error smoothing in all grid directions.

3. Matrix Structure: Structured n-Diagonal Matrices and Alternately Block Triangular Forms

Underlying the scheme is the exploitation of the sparsity and coupling reflected in the system matrix A arising from the discretization of the elliptic PDE on a structured grid (e.g., via a five-point stencil in 2D). A is characterized as a structured n-diagonal matrix, defined formally as:

$$\textbf{Definition (Structured n-diagonal matrix):}\quad \text{A square matrix } A \in \mathbb{R}^{n\times n} \text{ is called a structured n-diagonal matrix if its nonzero entries correspond exactly to a fixed computational molecule on a structured grid } Z \subset \mathbb{R}^d.$$

For the 2D case,

$A = G_1 + \beta I + G_2,$

where $G_1$ and $G_2$ are structured such that, under ordering induced by the sweeping schedule, one becomes block-lower-triangular and the other block-upper-triangular. These are instances of "alternatively block upper-lower triangular matrices," recursively defined via multilevel splittings. The significance of this structure is that the eigenvalues of G₁ and G₂ are simply the diagonal neighbors’ coefficients (e.g., $$a^{NE}_{i,j}, a^{NW}_{i,j}, a^{SE}_{i,j}, a^{SW}_{i,j}$$), which fundamentally simplifies convergence analysis.

4. Convergence Analysis and Relaxation Parameter Bounds

The convergence properties of the two-dimensional parallel approach are obtained through a block-wise splitting of A, leading to a system update across four directional sweeps per cycle. The iteration matrix T for a full cycle satisfies

$u^{(k+4)} = T u^{(k)} + (\text{RHS terms}),$

with an explicit bound

$$\|T\| \leq |1 - \omega_{SW}| \cdot |1 - \omega_{SE}| \cdot |1 - \omega_{NW}| \cdot |1 - \omega_{NE}|,$$

where each $\omega_{*}$ is the relaxation parameter for a given sweep direction. The convergence criterion reads

$$\prod_{d=1}^4 |1 - \omega_d| = |1-\omega_{SW}||1-\omega_{SE}||1-\omega_{NW}||1-\omega_{NE}| < 1,$$

and reduces to the classical SOR convergence bound $0 < \omega < 2$ if the relaxation parameters are uniform. This analysis generalizes familiar one-dimensional conditions to the multidirectional, decomposed parallel context, and is tightly linked to the block triangular structure of the underlying matrix.

Extensive numerical testing demonstrates that the parallelized method achieves iteration counts nearly identical to the sequential baseline, thus avoiding the convergence degradation seen in naïvely parallelized sweeps (Tavakoli, 2010).

5. Cache Efficiency, Data Locality, and Practical Implementation

A salient advantage of the sub-domain decomposition is enhanced data locality, which improves cache efficiency even in nominally sequential runs. By fitting sub-domains (and their associated data) within cache-resident memory, the algorithm sharply reduces cache-miss penalties relative to global sweeps. Each sub-domain prefetches necessary "halo" or boundary data, and then applies local updates with minimal external dependencies.

The efficiency factor is quantitatively defined as the ratio of cache-optimized run time to conventional sequential run time. Empirical results on large grids (e.g., with best performance for $2 \times 2$ or $3 \times 3$ sub-domain partitions in 2D/3D) validate both speedup and scalability, reinforcing the implementation viability.

The steps for practical implementation thus involve:

• Partitioning memory layout to support fast access for sub-domain sweeps
• Planning concurrency so that interface solves are synchronized without global barriers
• Adjusting sub-domain sizes for optimal cache utilization according to hardware specifics

6. Comparison with Alternative Parallelization Strategies

Conventional approaches to parallelizing sequential-sweep algorithms either fully decouple sub-domains (resulting in degraded convergence), or attempt synchronous updates with significant global communication overhead. In contrast, the overlapping-decomposition approach maintains local coupling via interface solves, which preserves both convergence properties and practical efficiency, even as the number of sub-domains increases.

In the language of domain decomposition methods, the presented method functions as a special-case overlapping Schwarz-type algorithm, but is distinctly adapted for structured grid sweeps and the associated matrix properties. The rigorous theoretical framework (matrix splitting, block-triangular structures, explicit convergence conditions) and empirical results together differentiate this two-dimensional parallel methodology from both classical and recent parallel smoothing/preconditioning schemes.

7. Broader Impact and Generalizations

The two-dimensional parallel approach pioneered in (Tavakoli, 2010) directly addresses the fundamental tension between maximizing parallel throughput and preserving algorithmic convergence in numerically sensitive, sweep-dependent solvers on structured grids. The framework is extensible to higher dimensions, irregular decompositions, and generalized n-diagonal matrix structures, provided the block splitting and interface-coupling mechanisms are appropriately adapted.

Real-world implications include:

• Enabling scalable, high-efficiency preconditioning for large elliptic PDE solvers in scientific and engineering simulations
• Supporting cache-efficient computations on modern architectures with deep memory hierarchies
• Providing a blueprint for parallelization of other inherently sequential sweeping paradigms, such as fast sweeping methods for Hamilton–Jacobi equations and certain iterative inversion schemes in computational imaging and geophysics

This innovative approach, grounded in both rigorous mathematics and careful systems-level optimization, offers a unified solution for high-performance computing tasks where sequential sweeps had previously been considered an inescapable bottleneck.