Parallel Banding Algorithm Plus (PBA+)

• The paper demonstrates that PBA+ extends traditional banded solvers by incorporating extraband nonzeros without compromising O(n) scaling.
• It employs a QP-decomposition with recursive forward and backward substitutions to achieve high numerical accuracy and efficient parallel implementation.
• Empirical results show that PBA+ outperforms standard Gaussian elimination in speed and error reduction for large-scale computational physics and chemistry problems.

The Parallel Banding Algorithm Plus (PBA+) is an $O(n)$ direct solver for large, sparse banded linear systems $A x = b$, where $A$ is an $n \times n$ matrix with nonzero entries confined within a band about the diagonal, and possibly with a small number $T_{\max}$ of extraband nonzero entries. PBA+ extends the core Parallel Banding Algorithm (PBA) framework by addressing these additional extraband elements without compromising linear scaling, and is designed for high accuracy and numerical efficiency, as required in computational physics and chemistry, such as finite-difference discretizations and holonomic constraint problems (García-Risueño et al., 2010).

1. Problem Definition and Matrix Structure

PBA+ targets systems where the matrix $A$ exhibits upper ($m_u$) and lower ($m_\ell$) half-bandwidths, such that

$$A_{i,j} = 0 \qquad \text{whenever}\ j - i > m_u \ \text{or}\ i - j > m_\ell.$$

Extra-band nonzero entries are incorporated as

$$A' = A + \sum_{t=1}^{T_{\max}} \delta_{i,R_t} \delta_{j,S_t} \Delta_t,$$

where each pair $A x = b$0 lies outside the original band, extending applicability to cases dominated by band structure with sparse, but non-negligible, extraband couplings (García-Risueño et al., 2010).

Key parameter notation includes:

• $A x = b$1, $A x = b$2,
• $A x = b$3,
• $A x = b$4.

2. Core Algorithmic Principle: QP-Decomposition and Recurrences

PBA+ exploits a unique decomposition, writing

$A x = b$5

where $A x = b$6 is lower-triangular with unit diagonal and $A x = b$7 is upper-triangular. The inverse is $A x = b$8. This structure leads to two substitution phases:

$A x = b$9

Constructive recursion builds $A$0, $A$1; each $A$2 operates on row $A$3, $A$4 on column $A$5, under explicit recursive updates for the coefficients $A$6 and $A$7:

• For $A$8 (diagonal and super-diagonal):

$A$9

$n \times n$0

• For $n \times n$1 (sub-diagonal):

$n \times n$2

In symmetric cases ($n \times n$3), the sub-diagonal recursion can be avoided by setting $n \times n$4.

After all $n \times n$5 are assembled, the solution proceeds via:

• Forward substitution for $n \times n$6:

$n \times n$7

• Backward substitution for $n \times n$8:

$n \times n$9

Each phase costs $T_{\max}$0, linear in $T_{\max}$1 for fixed bandwidths (García-Risueño et al., 2010).

3. Extension to Extra-Band Nonzeros: The Plus Variant

PBA+ generalizes the core banded approach by extending the recurrences for $T_{\max}$2 to all columns/rows occupied by extraband entries $T_{\max}$3. At each such location, the same analytic recurrences are applied—skipping zeros as appropriate—thus retaining linear scaling with a supplementary complexity term $T_{\max}$4, where $T_{\max}$5 is the number of extra-band entries and $T_{\max}$6 is a characteristic band width. This modular design facilitates accurate direct inversion for matrices that are banded plus a sparse pattern of extraband couplings without resorting to iterative refinement or full matrix storage (García-Risueño et al., 2010).

4. Parallelization Strategies

The core recursions couple only $T_{\max}$7 or $T_{\max}$8 preceding rows or columns, enabling efficient parallel execution via two main approaches:

a) Shared-memory (wavefront scheduling, Meurant '85 style):

• At each diagonal step $T_{\max}$9 ($A$0), all elements $A$1 (or $A$2) for $A$3 can be processed in parallel.
• Both forward and backward substitutions advance along anti-/diagonals, leading to a critical-path of $A$4, with parallel work $A$5 per step.
• On $A$6 cores: $A$7 (García-Risueño et al., 2010).

b) Distributed-memory (block-strip decomposition, Polizzi '06 style):

• The matrix is partitioned into $A$8 contiguous blocks, each handling local $A$9 evaluation up to block edges, and exchanging $m_u$0 boundary values.
• Each block performs a local factorization, then all blocks collectively solve a small coupling system of size $m_u$1.
• Final forward and backward solves require interface data broadcasts.
• Total time: $m_u$2; optimal for $m_u$3 (García-Risueño et al., 2010).

Near-linear speedup is achieved so long as $m_u$4 and the number of processors satisfies $m_u$5.

5. Computation of the Full Inverse

The factors $m_u$6 facilitate retrieval of the entire inverse $m_u$7 in $m_u$8 operations. Exploiting $m_u$9 and the banded nature, the entries are computed as:

$m_\ell$0

Each loop over the matrix involves $m_\ell$1 operations, resulting in $m_\ell$2 total complexity—competitive with dense inversion but leveraging full sparsity for memory efficiency (García-Risueño et al., 2010).

6. Complexity, Numerical Performance, and Comparative Table

PBA+ achieves direct $m_\ell$3 scaling for banded systems, maintaining leading constants that are empirically observed to be smaller than those for standard banded Gaussian elimination. For banded $m_\ell$4 sparse cases, the complexity is $m_\ell$5, with $m_\ell$6 as the number of extraband nonzeros. Storage requirements remain at $m_\ell$7. Forming the full inverse is $m_\ell$8 in both operations and storage.

Algorithm	FLOPs	Storage	Remarks
Gaussian elimination (banded)	$m_\ell$9	$$A_{i,j} = 0 \qquad \text{whenever}\ j - i > m_u \ \text{or}\ i - j > m_\ell.$$0	pivoting $$A_{i,j} = 0 \qquad \text{whenever}\ j - i > m_u \ \text{or}\ i - j > m_\ell.$$1 overhead
PBA	$$A_{i,j} = 0 \qquad \text{whenever}\ j - i > m_u \ \text{or}\ i - j > m_\ell.$$2	$$A_{i,j} = 0 \qquad \text{whenever}\ j - i > m_u \ \text{or}\ i - j > m_\ell.$$3	no LU storage; smaller constant
PBA+	$$A_{i,j} = 0 \qquad \text{whenever}\ j - i > m_u \ \text{or}\ i - j > m_\ell.$$4	$$A_{i,j} = 0 \qquad \text{whenever}\ j - i > m_u \ \text{or}\ i - j > m_\ell.$$5	$$A_{i,j} = 0 \qquad \text{whenever}\ j - i > m_u \ \text{or}\ i - j > m_\ell.$$6 = extra-band nonzeros
PBA full inverse	$$A_{i,j} = 0 \qquad \text{whenever}\ j - i > m_u \ \text{or}\ i - j > m_\ell.$$7	$$A_{i,j} = 0 \qquad \text{whenever}\ j - i > m_u \ \text{or}\ i - j > m_\ell.$$8

In extensive tests (up to $$A_{i,j} = 0 \qquad \text{whenever}\ j - i > m_u \ \text{or}\ i - j > m_\ell.$$9, $$A' = A + \sum_{t=1}^{T_{\max}} \delta_{i,R_t} \delta_{j,S_t} \Delta_t,$$0) on random banded matrices:

• With partial pivoting, PBA achieves $$A' = A + \sum_{t=1}^{T_{\max}} \delta_{i,R_t} \delta_{j,S_t} \Delta_t,$$1–$$A' = A + \sum_{t=1}^{T_{\max}} \delta_{i,R_t} \delta_{j,S_t} \Delta_t,$$2, nearly $$A' = A + \sum_{t=1}^{T_{\max}} \delta_{i,R_t} \delta_{j,S_t} \Delta_t,$$3-independent, with the error growing as $$A' = A + \sum_{t=1}^{T_{\max}} \delta_{i,R_t} \delta_{j,S_t} \Delta_t,$$4; Gaussian elimination yields $$A' = A + \sum_{t=1}^{T_{\max}} \delta_{i,R_t} \delta_{j,S_t} \Delta_t,$$5–$$A' = A + \sum_{t=1}^{T_{\max}} \delta_{i,R_t} \delta_{j,S_t} \Delta_t,$$6.
• PBA runtime is $$A' = A + \sum_{t=1}^{T_{\max}} \delta_{i,R_t} \delta_{j,S_t} \Delta_t,$$7 that of Gaussian elimination for $$A' = A + \sum_{t=1}^{T_{\max}} \delta_{i,R_t} \delta_{j,S_t} \Delta_t,$$8, and up to $$A' = A + \sum_{t=1}^{T_{\max}} \delta_{i,R_t} \delta_{j,S_t} \Delta_t,$$9 faster for $A x = b$00. Without pivoting, both algorithms deteriorate to $A x = b$01 error, but PBA remains $A x = b$02 faster.

For 1D Poisson problems ($A x = b$03, $A x = b$04), PBA completes in $A x = b$05 s with error $A x = b$06, compared to $A x = b$07 s for banded Gaussian elimination (Numerical Recipes), for a speedup of $A x = b$08.

In Lagrange multiplier enforcement for protein chains ($A x = b$09, $A x = b$10, $A x = b$11), PBA is $A x = b$12 faster than Gaussian elimination, and a symmetry-exploiting PBA+ variant further improves this to $A x = b$13, achieving errors down to $A x = b$14 (García-Risueño et al., 2010).

7. Implementation Considerations and Scalability

PBA+ relies on contiguous diagonal-oriented storage formats for bands, supporting cache-efficient sweeps of all $A x = b$15 elements needed. Extra-band entries are stored separately and invoked only in the extended recurrences. Partial pivoting, when adopted, incurs $A x = b$16 overhead and guards robustness against small $A x = b$17 elements.

Wavefront parallelism is effective for shared memory platforms as long as $A x = b$18; distributed memory strategies only require $A x = b$19 words exchanged at block interfaces, facilitating scalability for large $A x = b$20. In typical physical simulation scenarios, $A x = b$21 and $A x = b$22, so memory and computational costs remain highly favorable (García-Risueño et al., 2010).

Practical performance confirms PBA+'s suitability for large systems with predominantly banded structure and sparse extraband perturbations, offering robust accuracy, high efficiency, and broad parallel scalability.