Pith. sign in

REVIEW 3 major objections 5 minor 22 references

A Distributed-memory Tridiagonal Solver Based on a Specialised Data Structure Optimised for CPU and GPU Architectures

T0 review · 3 major / 5 minor · reviewed 2026-08-12 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read A distributed tridiagonal solver that needs only two neighbour-only MPI exchanges, and sustains near two-thirds of peak memory bandwidth on CPUs and GPUs.

desk verdict A solid performance-engineering paper whose distributed accuracy claim needs a proper error bound before the general 'diagonally dominant' promise is credible. read the letter →

arxiv 2411.13532 v1 pith:IZEIHV3Z submitted 2024-11-20 cs.DC physics.comp-ph

Semih Akkurt , Sébastien Lemaire , Paul Bartholomew , Sylvain Laizet This is my paper · ORCID
classification cs.DCphysics.comp-ph MSC 65F0565Y05
keywords tridiagonalsolverdistributedmemorydiagonaldominancecompactfinitedifferenceschemesMPIcommunicationGPUparallelisationbandwidthdatalayout
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper tries to show that solving batches of diagonally dominant tridiagonal systems on distributed-memory machines can be made both communication-light and bandwidth-efficient. It introduces DistD2-TDS, an algorithm whose communication cost does not grow with the number of ranks: after a local decoupling phase, the global system reduces to independent $2\times2$ systems across each MPI boundary, so the solver only exchanges halo data and one set of boundary values with immediate neighbours. A companion data structure packs a small group of tridiagonal lines contiguously in memory, which makes the sequential passes vectorisable on CPUs and thread-parallel on GPUs and keeps data movement close to the practical memory-bandwidth ceiling. If the claims hold, compact-scheme PDE solvers that currently pay for global transposes or multi-step reductions can instead scale to thousands of ranks with neighbour-only communication while sustaining around two-thirds of theoretical peak bandwidth.

What carries the argument

The specialised data structure packs tridiagonal lines in groups of size $SZ$, typically 8 for CPU vector registers and 32 for GPU streaming multiprocessors, with the $n$-th point of each system adjacent in memory; this layout is what makes vectorisation, thread-level parallelism, cache blocking, and fused right-hand-side construction possible. The algorithmic mechanism is the parallel-diagonally-dominant reduction: local inversion decouples each subdomain's interior from its boundaries, and the resulting penta-diagonal coupling is truncated to independent $2\times2$ systems across MPI boundaries because diagonally dominant systems make the coupling terms fall below the zero-machine value.

What would settle it

Run the paper's sixth-order compact first-derivative solve, with the implicit coefficient $\alpha$ in Equation 2 pushed toward $1/2$ so the system is only weakly diagonally dominant, using subdomains of a few points per rank across many ranks; if the maximum error relative to the serial solver grows as dominance weakens, subdomains shrink, or ranks increase, the discarded coupling terms are not negligible.

Watch

Extended reading notes

Core claim

DistD2-TDS solves a batch of diagonally dominant tridiagonal systems in three phases—decoupling, neighbour communication, substitution—and the decoupling phase transforms the global problem so that only $2\times2$ systems across subdomain boundaries remain. The central claim is that the reduced penta-diagonal system produced by multiplying each local region by its inverse contains entries below the zero-machine value for diagonally dominant systems; discarding them leaves independent $2\times2$ systems, so solving the whole distributed problem requires only two sets of MPI communications, a halo exchange and then a boundary exchange, with the number of exchanges independent of rank count. The accompanying data structure stores $SZ$ tridiagonal systems interleaved point-by-point, giving contiguous, predictable memory access and enabling cache blocking and kernel fusion; benchmarks show the solver's runtime tracks STREAM-style copy/scale data-movement costs, with 66% to 68% of theoretical peak bandwidth sustained on a realistic three-dimensional PDE at scale.

Load-bearing premise

Everything rests on the claim that for a diagonally dominant tridiagonal system the coupling terms left after the local decoupling are smaller than the floating-point zero and can be thrown away without changing the answer; the paper's accuracy test covers only one compact first-derivative scheme with 32 to 256 points on two ranks.

Editorial extensions

If this is right

  • Communication per distributed solve becomes independent of rank count: only a halo exchange and one nearest-neighbour boundary exchange are needed, so strong scaling does not degrade as ranks are added.
  • For compact-scheme PDE solvers, the per-step cost of evaluating derivatives approaches the cost of reading and writing the field once, near the practical memory-bandwidth ceiling—66% to 68% of theoretical peak in the paper's large-scale runs.
  • The layout-reordering overhead of switching between spatial directions in a 3D PDE is small, at 7% to 12% of total data movement, so the strategy is viable for full 3D solvers rather than only isolated line solves.
  • The same diagonal-dominance reduction can in principle be extended to generic non-diagonally-dominant tridiagonal systems, as the paper states, which would widen the class of solvers that avoid global communication.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • If the communication bound is real, then for a fixed subdomain size the solver's runtime should be nearly flat as ranks increase, with only halo depth and memory bandwidth left as scaling costs; a weak-scaling sweep at fixed points-per-rank would test this directly.
  • The accuracy of the discarded coupling terms is demonstrated for only one compact scheme with 32 to 256 points on two ranks; a systematic sweep over diagonal-dominance strength, subdomain size, and rank count would show where the 'below zero-machine value' assumption starts to fail and could yield a quantitative safety criterion for DistD2-TDS.
  • The same packed-group data layout could be reused for other bandwidth-bound stencil or banded solvers, such as pentadiagonal or block-tridiagonal systems, potentially extending the vectorisation and thread-level-parallelism benefits beyond tridiagonal solves.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 5 minor

Summary. The paper proposes DistD2-TDS, a distributed-memory tridiagonal solver for batches of diagonally dominant tridiagonal systems. The algorithm combines a specialized data layout that packs groups of tridiagonal systems contiguously in memory with a distributed strategy based on the Parallel Diagonal Dominant (PDD) idea: after local decoupling phases, the reduced global penta-diagonal system is claimed to collapse into independent 2x2 systems across MPI boundaries, requiring only two sets of neighbor communications. The paper reports single-rank performance on CPUs and GPUs, strong scaling up to 384 NVIDIA H100 GPUs and 8192 AMD EPYC 7742 CPUs, and an application to a 3D non-linear PDE, with sustained bandwidth around 66% of theoretical peak. The central accuracy claim rests on the assertion that certain entries of the reduced system are 'below the zero-machine value' for diagonally dominant systems, inherited from Sun (1995), but no error bound or parametric study is provided for DistD2.

Significance. If the correctness claim holds, the paper offers a significant advance: communication cost that does not grow with the number of ranks (two neighbor exchanges per solve), combined with a data structure that achieves near-peak bandwidth on modern CPUs and GPUs. The performance analysis is carefully anchored to BabelStream copy/scale benchmarks, with explicit data-movement accounting in Tables 3-5, and the large-scale benchmarks on up to 8192 CPU cores and 384 H100 GPUs are valuable. The practical 3D PDE demonstration, including kernel fusion and reordering overheads, is a strength. However, the correctness of the distributed algorithm rests on an unquantified approximation whose validity range is not established; the current accuracy evidence is limited to a single compact-scheme family at small sizes and two ranks. The paper's advertised scope of general diagonally dominant tridiagonal systems is therefore not yet supported.

major comments (3)
  1. [Section 2.2.2 and Section 4.2, Figure 8, Figure 10] The central algorithmic claim that the reduced penta-diagonal system collapses to independent 2x2 systems relies on the assertion that certain entries are 'below the zero-machine value' for diagonally dominant systems, inherited from Sun [7]. No bound is derived in this manuscript for DistD2 quantifying the error as a function of the diagonal-dominance ratio (rho) and subdomain length (m). The only accuracy test (Figure 10) covers a single sixth-order compact first-derivative scheme with alpha=1/3 (dominance ratio 2/3), n=32-256, and 2 ranks; it does not sweep dominance strength, subdomain size, or rank count. For weakly dominant systems or short subdomains, the dropped couplings decay only like rho^m and can easily exceed machine epsilon, causing silent accuracy loss. The abstract and introduction generalize to 'diagonally dominant tridiagonal systems,' so the paper needs either a quantitative error bound or a condition relating rho and m, together with validation across that parameter range.
  2. [Algorithms 5-7] The pseudocode is not sufficiently specified for verification. In Algorithm 5, line 18 computes w1 = 1/(1 - sc_1 sa_2 with a missing closing parenthesis; r_j is defined only for j>=3 in the loop at lines 7-12, yet Algorithm 6 uses r_j for j=1:2 (r1 and r2 are undefined). Algorithm 7 uses d_i^s and d_i^e without defining where these values come from (presumably from the 2x2-system communication phase, but the connection is not stated). The substitution loop in Algorithm 7 iterates over j=1:n and then overwrites the boundary entries di,1 and di,n; this is correct only if the boundary entries are not meant to be updated by the loop, so the loop should be restricted to j=2:n-1 or the overwrite should be stated explicitly. These issues impede reproduction and should be fixed.
  3. [Section 1 and Section 4] The introduction states that 'limitations of our strategy and potential extensions are discussed in detail in Section 4,' but Section 4 (Proposed Strategy) contains no such discussion. The paper should either include a limitations subsection in Section 4 or correct the pointer. More substantively, the limitations of the diagonal-dominance assumption - including the valid range of dominance strength, the effect of subdomain size, and the behavior for nearly non-dominant systems - should be discussed explicitly, as they are directly relevant to the central correctness claim.
minor comments (5)
  1. [Section 5.2] The strong-scaling efficiencies (82.7%, 95.8%, etc.) are reported without defining the baseline or the formula used; specify whether they are relative to the runtime at the smallest node count and whether the problem size is held constant.
  2. [Algorithm 5, lines 1-6] The preprocessing loop for j=1:2 assumes input arrays sa_j and sc_j are already initialized, but this is not stated; clarify that these are the input matrix coefficients.
  3. [Algorithm 7, lines 4-8] The loop 'for j = 1 : n' is later followed by assignment to di,1 and di,n; consider changing the loop bound to 'j = 2 : n-1' to avoid the appearance of overwriting boundary values.
  4. [Figure 6 caption] The caption says 'Data continuity in memory is in column-major order,' but the figure itself is a schematic of the grouping; clarify the notation (SZ, nx, ny*nz/SZ) and how it maps to the displayed domain.
  5. [Table 4] The peak bandwidth values for AMD EPYC 7742 are listed as 205 GB/s, while Table 1 reports a newer EPYC 9754 with 460.8 GB/s; ensure the architectural context is clear so readers do not confuse the two.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: DistD2-TDS is not defined in terms of the quantities it predicts, and its correctness and performance claims are anchored to external references (Thomas solve and BabelStream).

full rationale

The paper's central algorithmic claim is that a diagonally dominant tridiagonal system can be reduced to independent 2x2 systems across MPI boundaries, requiring only neighbour communication. This claim is not circular: it is derived from the structure of the decoupling/substitution phases in Algorithms 5-7 and is justified by an external reference, Sun (1995), not by a self-citation. The accuracy of dropping 'below zero-machine' entries is then checked directly against the Thomas algorithm in Figure 10, so the claimed result is not defined in terms of itself. The accuracy test is narrow (one compact scheme, n=32-256, 2 ranks), and the Introduction promises a limitations discussion in Section 4 that is not actually present; these are validation and documentation gaps, not circular reasoning. The performance claims are also non-circular: bandwidth utilisation in Table 6 is computed from measured runtimes divided by independently enumerated data-movement requirements, with BabelStream as an external bandwidth baseline. The parameter SZ is chosen from hardware vector/thread widths, not fitted to force the reported bandwidth numbers. The self-citations present (Incompact3d as a use-case example, a cache-blocking reference for data packing) are contextual and not load-bearing. No fitted parameter is renamed as a prediction, and no equation is equivalent to its own input by construction. The headline results therefore have independent content and the proper finding is no significant circularity, score 0.

Assumptions & free parameters 1 free parameters · 4 assumptions · 0 invented entities

The correctness of DistD2 rests on the PDD approximation theorem (axiom 2), which is only cited and not bounded by the authors. The performance claims rest on the hand-tuned block size SZ and on the BabelStream ceiling assumption. No new physical entities are introduced; the packed data structure is the paper's contribution, not an unverified postulate.

free parameters (1)
  • SZ (group size) = 8 (CPU), 32 (GPU)
    Number of tridiagonal systems packed together in the specialised data structure. Chosen from vector register width (512-bit AVX) and SM double-precision core count; not swept in the paper, yet it directly affects the measured bandwidth and runtime.
assumptions (4)
  • standard math Thomas algorithm solves tridiagonal systems exactly (Algorithm 1).
    Used as the reference method in the accuracy test (Figure 10) and as the starting point for the decoupling phase in Section 4.2.
  • domain assumption For diagonally dominant tridiagonal matrices, the coupling terms outside each 2x2 boundary block in the reduced system are below machine zero and can be dropped.
    Imported from Sun 1995 (PDD); the paper does not state the quantitative conditions on diagonal dominance or subdomain size, and does not prove them for the modified-Thomas variant. Invoked in Section 2.2.2 and Figure 8.
  • domain assumption Compact schemes of the form of Equation (2) with coefficients from Lele (1991) produce strictly diagonally dominant tridiagonal systems.
    The paper states this in Section 1 and relies on it for the application; for alpha = 1/3 it holds by inspection since |a| + |c| = 2/3 < 1 = |b|, but the general statement over various numerical methods is asserted, not proven.
  • domain assumption BabelStream copy and scale benchmarks represent the practical attainable bandwidth ceiling of the hardware (75-85% of peak).
    The performance model in Section 5.1 treats these benchmarks as the reference ceiling; the claim is taken from Deakin et al. (BabelStream), not verified in this paper.

how reviews work

0 comments
Cite this review

Pith. "Pith review of A Distributed-memory Tridiagonal Solver Based on a Specialised Data Structure Optimised for CPU and GPU Architectures." pith.science (2026). https://pith.science/paper/IZEIHV3Z

@misc{pith2026241113532,
  author       = {Pith},
  title        = {Pith review of: A Distributed-memory Tridiagonal Solver Based on a Specialised Data Structure Optimised for CPU and GPU Architectures},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/IZEIHV3Z}},
  note         = {Machine review of arXiv:2411.13532}
}
read the original abstract

Various numerical methods used for solving partial differential equations (PDE) result in tridiagonal systems. Solving tridiagonal systems on distributed-memory environments is not straightforward, and often requires significant amount of communication. In this article, we present a novel distributed-memory tridiagonal solver algorithm, DistD2-TDS, based on a specialised data structure. DistD2-TDS algorithm takes advantage of the diagonal dominance in tridiagonal systems to reduce the communications in distributed-memory environments. The underlying data structure plays a crucial role for the performance of the algorithm. First, the data structure improves data localities and makes it possible to minimise data movements via cache blocking and kernel fusion strategies. Second, data continuity enables a contiguous data access pattern and results in efficient utilisation of the available memory bandwidth. Finally, the data layout supports vectorisation on CPUs and thread level parallelisation on GPUs for improved performance. In order to demonstrate the robustness of the algorithm, we implemented and benchmarked the algorithm on CPUs and GPUs. We investigated the single rank performance and compared against existing algorithms. Furthermore, we analysed the strong scaling of the implementation up to 384 NVIDIA H100 GPUs and up to 8192 AMD EPYC 7742 CPUs. Finally, we demonstrated a practical use case of the algorithm by using compact finite difference schemes to solve a 3D non-linear PDE. The results demonstrate that DistD2 algorithm can sustain around 66% of the theoretical peak bandwidth at scale on CPU and GPU based supercomputers.

Figures

Figures reproduced from arXiv: 2411.13532 by the authors.

Figure 1
Figure 1. Thomas algorithm for solving tridiagonal matrice [PITH_FULL_IMAGE:figures/full_fig_p006_1.png] view at source ↗
Figure 2
Figure 2. 1D/2D/3D decomposition strategies with all their [PITH_FULL_IMAGE:figures/full_fig_p008_2.png] view at source ↗
Figure 3
Figure 3. This allows distributed-memory algorithms to hav [PITH_FULL_IMAGE:figures/full_fig_p008_3.png] view at source ↗
Figures from the paper (11 more)
Figure 3
Figure 3. Figure 3: Subdomains in a distributed-memory tridiagonal m [PITH_FULL_IMAGE:figures/full_fig_p009_3.png]
Figure 4
Figure 4. Figure 4: First phase of the modified Thomas algorithm. Each s [PITH_FULL_IMAGE:figures/full_fig_p010_4.png]
Figure 5
Figure 5. Figure 5: Parallel diagonally dominant algorithm. Each region is multiplied by its local inverse which results in a reduced penta-diagonal system with first and last data points from each rank. 12 [PITH_FULL_IMAGE:figures/full_fig_p012_5.png]
Figure 6
Figure 6. Figure 6: The proposed data structure for an x-directional tridiagonal system. Data continuity in memory is in column-major order. FLOPs per cycle, and for a double precision simulation on GPUs, SZ = 32 as a single streaming multiprocessor (SM) typically has 32 double precision …
Figure 7
Figure 7. Figure 7: A single group with SZ = 4 and n = 15. The tridiagonal system is demonstrated with a single RHS vector, and the location of the RHS vector in the specialised data structure is highlighted in colour to indicate the exact distribution of a single RHS in memory. The data …
Figure 8
Figure 8. Figure 8: Distributed-memory tridiagonal matrix algorith [PITH_FULL_IMAGE:figures/full_fig_p024_8.png]
Figure 9
Figure 9. Figure 9: Visualisation of the fused forward pass and RHS con [PITH_FULL_IMAGE:figures/full_fig_p024_9.png]
Figure 10
Figure 10. Figure 10: Order of accuracy of DistD2 algorithm compared ag [PITH_FULL_IMAGE:figures/full_fig_p026_10.png]
Figure 11
Figure 11. Figure 11: Single rank performance of the DistD2 and Thomas a [PITH_FULL_IMAGE:figures/full_fig_p029_11.png]
Figure 12
Figure 12. Figure 12: Single rank performance of the DistD2 and Thomas a [PITH_FULL_IMAGE:figures/full_fig_p029_12.png]
Figure 13
Figure 13. Figure 13: Strong scaling of DistD2 algorithm on ARCHER2 [19 [PITH_FULL_IMAGE:figures/full_fig_p034_13.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

22 extracted references · 22 canonical work pages

  1. [7]

    Sun, Application and accuracy of the parallel diag onal dominant algorithm, Parallel computing 21 (8) (1995) 1241–1267

    X.-H. Sun, Application and accuracy of the parallel diag onal dominant algorithm, Parallel computing 21 (8) (1995) 1241–1267

    work page 1995

  2. [1]

    S. K. Lele, Compact finite difference schemes with spectra l-like resolu- tion, Journal of computational physics 103 (1) (1992) 16–42

    work page 1992

  3. [2]

    D. W. Peaceman, H. H. Rachford, The numerical solution of parabolic and elliptic differential equations, Journal of the Society for Industrial and Applied Mathematics 3 (1) (1955) 28–41

    work page 1955

  4. [3]

    László, M

    E. László, M. Giles, J. Appleyard, Manycore algorithms f or batch scalar and block tridiagonal solvers, ACM transactions on mathema tical soft- ware 42 (4) (2016) 1–36

    work page 2016

  5. [4]

    G. D. Balogh, T. S. Flynn, S. Laizet, G. R. Mudalige, I. Z. R eguly, Scal- able many-core algorithms for tridiagonal solvers, Comput ing in science & engineering 24 (1) (2022) 26–35

    work page 2022

  6. [5]

    Kim, J.-H

    K.-H. Kim, J.-H. Kang, X. Pan, J.-I. Choi, PaScaL_TDMA: A library of parallel and scalable solvers for massive tridiagonal syst ems, Computer physics communications 260 (2021) 107722–. 40

    work page 2021

  7. [6]

    H. Song, K. V. Matsuno, J. R. West, A. Subramaniam, A. S. Gh ate, S. K. Lele, Scalable parallel linear solver for compact band ed systems on heterogeneous architectures, Journal of Computational Physics 468 (2022) 111443

    work page 2022

  8. [8]

    L. H. Thomas, Elliptic problems in linear difference equa tions over a network, Watson Sc. Comp. Lab. Rep., Columbia University, N ew York (1949)

    work page 1949

Show all 22 references
  1. [9]

    Sherman, W

    J. Sherman, W. J. Morrison, Adjustment of an Inverse Matr ix Corre- sponding to a Change in One Element of a Given Matrix, The Anna ls of Mathematical Statistics 21 (1) (1950) 124 – 127

    1950

  2. [10]

    Laizet, N

    S. Laizet, N. Li, Incompact3d: A powerful tool to tackle turbulence problems with up to O(105) computational cores, International Journal for Numerical Methods in Fluids 67 (11) (2011) 1735–1757

    2011

  3. [11]

    Rolfo, C

    S. Rolfo, C. Flageul, P. Bartholomew, F. Spiga, S. Laize t, The 2DE- COMP&FFT library: an update with new CPU/GPU capabilities, Jour- nal of Open Source Software 8 (91) (2023) 5813

    2023

  4. [12]

    Gander, G

    W. Gander, G. H. Golub, Cyclic reduction - history and ap plications, in: In Proceedings of the Workshop on Scientific Computing, 1 997

  5. [13]

    Akkurt, F

    S. Akkurt, F. Witherden, P. Vincent, Cache blocking str ategies applied to flux reconstruction, Computer Physics Communications 27 1 (2022) 108193. 41

    2022

  6. [14]

    OpenMP Reference Guide, https://www.openmp.org/res ources/refguides/ (2021)

    2021

  7. [15]

    CUDA Programming Model, https://docs.nvidia.com/cu da/cuda-c- programming-guide/ (2024)

    2024

  8. [16]

    J. D. McCalpin, Memory bandwidth and machine balance in current high performance computers, in: IEEE Computer Society Tech nical Committee on Computer Architecture (TCCA) Newsletter, 199 5, p. 19–25

  9. [17]

    J. D. McCalpin, Stream2, https://www.cs.virginia.ed u/stream/stream2/ (1999)

    1999

  10. [18]

    Deakin, J

    T. Deakin, J. Price, M. Martineau, S. McIntosh-Smith, E valuating at- tainable memory bandwidth of parallel programming models v ia babel- stream, International Journal of Computational Science an d Engineering 17 (3) (2018) 247–262

    2018

  11. [19]

    ARCHER2 Hardware Specifications, https://www.archer 2.ac.uk/about/hardware.html (2021)

    2021

  12. [20]

    ac.uk/about/hardware.html

    Cirrus Hardware Specifications, https://www.cirrus. ac.uk/about/hardware.html

  13. [21]

    Isambard 3 Hardware Specifications, https://docs.isambard.ac.uk/specs/#system-specificat ions-isambard-3- grace (2024)

    2024

  14. [22]

    Marenostrum 5 Hardware Specifications, https://www.bsc.es/ca/marenostrum/marenostrum-5 (202 3). 42

Pith tools

Reviewed August 12, 2026 · model on record in the stance chip above.