ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE RAPHAEL S. STEINER, CHRISTOS K. MATZOROS, PÁL ANDRÁS PAPP, TONI BÖHNLEIN, AND ALBERT-JAN N. YZELMAN
arXiv:2607.02324v1 [cs.DC] 2 Jul 2026
Abstract. We introduce stale synchronous parallel as a mode of execution in parallel sparse triangular linear system solve and present a general directed-acyclic-graph scheduler capable of producing such schedules. Stale-synchronous-parallel schedules allow the overlap of synchronisation and compute which results in a geometric-mean speed-up of 7-30% of our scheduler, ElasticDivide, over state-of-the-art synchronous scheduler GrowLocal on an ARM machine using 48 cores. On an x86 machine using 48 cores, we report geometric-mean speedups of 19-60% over SpMP.
1. Introduction Sparse triangular solve is an omnipresent operation in computing. Whether it be in engineering, data analytics, artificial intelligence, or scientific computing, systems of equations need to be solved. This typically involves (sparse) triangular solve as part of the solving process directly, through methods like LU, QR, and Cholesky factorisations and Gauß–Seidel, or indirectly through (pre-)conditioning of the system for faster iterative solving methods such as conjugate gradient. When solving ever bigger problem instances, sparsity plays an evermore important role. It allows one to cut redundant compute, but this comes at the cost of losing structure. This results in fine-grained dependencies and operations which makes sparse triangular solve (SpTrSV) a hard problem to parallelise. Therefore, getting the most performance out of modern multi-core architectures is a challenging endeavour. To combat this, SpTrSV performance is often optimised in larger contexts such as linear or symmetric solves where preprocessing allows for the reintroduction of structure. A prominent example of this is the nested dissection technique [Geo73, LRT79, KK98, APc04, GBDD10, BAvL+ 19] which concentrates non-zeroes for locality and introduces coarse-grained parallelism [DS05, LDS+ 23]. Another important technique is the introduction of so-called supernodes [CAGL+ 87, Li05, FZW+ 23, LNP93, YRE20, SG04, SGFS01]. Supernodes give rise to several optimisations. First, the grouping of operations into supernodes allows for the use of highly optimised dense kernels, and second, the dependency graph becomes coarser which enables dynamically computing a parallel schedule or dynamic dispatching methods as their benefits now outweigh their overhead. Another notable and generally applicable technique is the (recursive) splitting of SpTrSV into two SpTrSV problems of half the size and an easier-to-parallelise SpMV in-between [AS89, May09, LLH+ 16, LNL20, AYU21]. After applying these structural techniques (if they are feasible), one is left with a parallel scheduling problem. The tasks to be scheduled typically have irregular dependencies and their size can range from whole blocks or supernodes to a couple of floating-point operations. It is this scheduling problem that we address in this paper. More precisely, we consider the fine-grained scheduling problem of parallelising the forward-substitution algorithm, cf. Algorithm 2.1. Nevertheless, our algorithm and the involved ideas may, of course, be applied to any parallel scheduling problem. To get the most performance out of SpTrSV, a parallel scheduling algorithm must • balance workloads among cores, • limit coordination overhead, and • consider spatial and temporal locality. Key words and phrases. Sparse triangular linear system solve, SpTrSV, SpTrSM, forward- and backwardsubstitution algorithm, stale-synchronous-parallel algorithm. 1
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE HDagg
GrowLocal
ElasticDivide
AMD iCholesky
Metis ND
40
ARM Kunpeng
Speed-up over Serial
SuiteSparse
2
30 20 10 0
16
32
48
64
Cores
80
96
16
32
48
64
Cores
80
96
16
32
48
64
80
96
Cores
Figure 1.1. Geometric-mean speed-ups over Serial for various data sets, cf. §5.2, on ARM Kunpeng with 10-th to 90-th percentile range shown. Several algorithms to this end have been proposed in the literature. Broadly speaking, they come in two flavours, synchronous and asynchronous, depending on the type of schedule they produce. Early algorithms include asynchronous self-scheduling [SMB88] and the synchronous so-called wavefront schedulers [AS89, Sal90]. Both approaches suffered from a lot of synchronisation overhead due to frequent synchronisation [RG92, PSSD14]. In a breakthrough paper [PSSD14], Park et al. addressed this issue. Their scheduler SpMP combines the coarseness from the grouping of wavefront schedulers with asynchronous compute and point-to-point synchronisation. To further reduce the number of synchronisations, they employed an approximate transitive reduction to get rid of the majority of implied (and thus unnecessary) dependencies. Yılmaz et al. [YSAU20] were further able to reduce the number of dependencies by forcing the cores to be out-of-sync by only a bounded amount. Additionally, their scheduler ALB takes into account the amount of parallelism available locally in the matrix as the computation moves forward to decide how many of the available cores to use. On the synchronous front, there have also been improvements. The HDagg scheduler by Zarebavani et al. [ZCL+ 22] combines successive wavefronts whenever beneficial in order to reduce the number of synchronisations. In contrast, Böhnlein et al. [BPSY25, BPS+ 25] moved away from wavefronts and took inspiration from list scheduling [Gra69, ACD74, HCAL89, RVG02, MSQ03]. Their scheduler, GrowLocal, interleaves local list scheduling with synchronisation barriers which communicate the completion of rows. This allowed for much longer independent compute, reducing the number of synchronisations by an order of magnitude. 1.1. Our contribution. In this paper, we strike a middle ground between synchronous and asynchronous compute. Namely, we employ for the first time stale synchronous parallel [CHK+ 13, CCH+ 14] as the mode of execution in sparse triangular solve. This execution model explicitly overlaps synchronisation and compute by forcing inter-core dependencies to be elongated. In other words, the information computed by one core can only be used on another core after several synchronisation events. Since the information is not immediately required after a synchronisation event, the synchronisation can be overlapped with other compute based on already locally available information. In the case of synchronisation-heavy workloads with many threads, we show that this can result in a 2× geometric-mean speed-up over synchronous execution. On architectures with higher number of NUMA-domains, this can go up to 4.5×. Even when the architecture is more uniform and the synchronisation events are few, it leads to improvement in the double digit percentages. Naturally, stale synchronous parallel puts additional strain on the schedule to be generated and therefore also on the scheduler. Nevertheless, we demonstrate empirically that the computational graphs associated to an SpTrSV of a matrix are elastic enough to allow for such schedules. Indeed, our scheduler ElasticDivide produces about the same number of synchronisation barriers (0.85× to 1.28×) as the state-of-the-art synchronous scheduler GrowLocal [BPS+ 25] when appropriately normalised, yet overlaps compute with synchronisation. In Figure 1.1, we depict geometric-mean speed-ups over Serial on an ARM Kunpeng machine with 48 cores and two NUMA-domains per socket. The three plots correspond to three data sets which are derived from the SuiteSparse Matrix Collection [DH11]. SuiteSparse consists
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
3
of matrices directly taken from the data set, and the other two data sets are application oriented with the matrices having been preprocessed with an approximate-minimal-degree ordering and incomplete Cholesky using Eigen [GJ+ 10], respectively with nested dissection using Metis [KK98]. The reported geometric-mean speed-ups range from 3% on small core counts to over 100% on large core counts over the second best measured algorithm GrowLocal. On 48 cores (single socket), the geometric-mean speed-ups are 30%, 7%, and 17%, respectively. SpMP, being x86-specific, is missing from Figure 1.1. On 48 cores on an AMD EPYC machine, we achieve geometric-mean speed-ups of 19%, 60%, and 26%, respectively, over SpMP with geometric-mean speed-ups reaching over 100% on smaller core counts. In summary, in this paper, we • propose and analyse solving SpTrSV problems via stale synchronous parallelisation, • introduce ElasticDivide, a stale-synchronous-parallel scheduler for general directed acyclic graphs, and • show their effectiveness on non-uniform memory architectures. 1.2. Overview. We describe the stale-synchronous-parallel scheduling problem in §2. The ElasticDivide scheduler is presented in detail in §3. In §4, we describe the implementation of the kernel and the barrier. In §5 and §6, we describe our evaluation procedure and evaluate the algorithms, followed by a discussion in §7. Acknowledgements. We would like to thank Weifeng Liu, Olaf Schenk, Xiaoye Li, Dimosthenis Pasadakis, Lorenzo Migliari, and Kiril Dichev for stimulating conversations on this and surrounding topics. 2. Preliminaries 2.1. Problem description. In sparse triangular solve (vector), one is given a sparse invertible matrix A = (Ai,j )i,j=1,...,n ∈ Rn×n and a vector b = (bi )i=1,...,n ∈ Rn and seeks to compute the vector x = (xi )i=1,...,n ∈ Rn such that Ax = b
(2.1)
holds. The standard algorithm to solve for x is the forward-substitution algorithm if the matrix A is lower triangular, respectively the backward-substitution algorithm if A is upper triangular. Assume from now on that A is lower triangular. The forward-substitution algorithm is described in Algorithm 2.1. Algorithm 2.1: Forward-substitution algorithm Data: An invertible lower-triangular matrix A ∈ Rn×n and a vector b ∈ Rn . Result: A vector x ∈ Rn such that Ax = b. 1 for i = 1, . . . , n do 2 t ← bi 3 4
for j = 1, . . . , i − 1 with Ai,j ̸= 0 do t ← t − Ai,j xj
5
xi ← t/Ai,i
In contrast to the dense case, when the matrix A is sparse, there is opportunity to parallelise the outer for-loop, cf. Line 1. Though, one has to be careful as there are dependencies. The i-th iteration depends directly on the j-th iteration if and only if Ai,j ̸= 0 and i ̸= j. More generally, the i-th iteration depends on the j-th iteration if and only if there is a sequence i = ℓ0 > ℓ1 > · · · > ℓm = j with Aℓk ,ℓk+1 ̸= 0 for k = 0, . . . , m − 1 for some m ≥ 1. We can and shall capture this in a directed acyclic graph GA = (VA , EA ), where VA = {1, . . . , n}, and
(2.2)
EA = {(j, i) ∈ VA × VA | Ai,j ̸= 0} .
(2.3)
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
4
We may also attach a (vertex-)weight to GA modelling the execution time of the i-th iteration corresponding to a vertex i ∈ VA . A natural choice as this weight is the number of non-zero elements the corresponding row, that is ωA : VA → Z>0 i 7→ ♯{j ∈ VA | Ai,j ̸= 0} .
(2.4)
This is also equal to the in-degree plus one in the graph representation. From here on out, we work with the graph representation G = GA of the sparse lower triangular matrix A and formulate the problem as a scheduling problem in an execution model described in the next section. 2.2. Stale synchronous parallel. Stale synchronous parallel is a computational model extending the bulk-synchronous-parallel model [Val90a, Val90b]. It was formally introduced by Cui et al. [CCH+ 14] and originally stems from a paper by Cipar et al. [CHK+ 13] in machine learning. Similar to the bulk synchronous parallel, stale synchronous parallel operates in ‘rounds’ which are called supersteps, but, unlike bulk synchronous parallel, cores may move the computation ahead into subsequent supersteps as long as the maximal difference in the superstep index remains bounded by a fixed quantity called the staleness. Formally, a compute schedule in the stale-synchronous-parallel model may be defined as follows. Definition 2.1. A stale-synchronous-parallel schedule of staleness s ∈ Z>0 of a directed acyclic graph G = (V, E) to a set of cores P consists of a pair (π, σ) of maps π : V → P and σ : V → Z≥0 satisfying ∀(v, w) ∈ E : σ(v) + s · δπ(v),π(w) ≤ σ(w) ,
(2.5)
where δi,j is the Kronecker delta. The map π maps vertices to cores and the map σ maps vertices to supersteps. Remark 2.2. The definition of a stale-synchronous-parallel schedule of staleness s = 1 is the same as a bulk-synchronous-parallel schedule. Remark 2.3. A stale-synchronous-parallel schedule of staleness s is naturally a stale-synchronous-parallel schedule of staleness t such that 1 ≤ t ≤ s, in particular also a bulk-synchronousparallel schedule. The use of the word ‘stale’ in the name ‘stale synchronous parallel’ is perhaps unfortunate and stems from the first paper [CHK+ 13], where it was used in machine-learning model training. Workers were allowed to compute model updates based on a mildly outdated (stale) model, which allowed for some asynchronicity whilst (provably) maintaining convergence and quality guarantees. Only later [CCH+ 14], it was realised that the ‘data’ need not be outdated (or stale), but instead having ‘long’ dependencies suffices. The stale-synchronous-parallel model has since found many applications, such as recommendation systems [DPY25], matrix factorisation [HCC+ 13], topic modeling [HCC+ 13], and training for deep learning [ZALC19]. We also note that overlapping computation and communication, a key motivation for the stale synchronous parallel model, had already been identified as a core feature of the bulk synchronous parallel model [Val90a]. 3. The ElasticDivide scheduler Our scheduling algorithm ElasticDivide is an extension of the algorithm GrowLocal presented in [BPS+ 25]. A sketch of ElasticDivide is given by Algorithm 3.1. The main novelty of our algorithm is the adaptation to a stale-synchronous-parallel schedule of staleness 2, cf. Definition 2.1. To achieve this, we had to innovate and add two key features: (i) a delay in resolving dependencies across cores, and (ii) strike a balance between assigning vertices in the current superstep and the next.
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
5
Algorithm 3.1: Sketch of ElasticDivide scheduler Data: A vertex-weighted directed acyclic graph G = (V, E, ω) and cores P = {1, 2, . . . , k}. Result: Assignments of vertices to cores π : V → P and superstep σ : V → Z≥0 . Rule I: A vertex v is assignable to core p and superstep s ∈ Z≥0 if and only if ∀(w, v) ∈ E : σ(w) + 2 · δπ(w),p ≤ s . Rule II: Vertices are prioritised according to (i) core exclusivity, and then (ii) smallest ID. 1 s←0 2 while not all vertices are assigned do 3 α ← 20 4 5 6 7 8 9 10 11 12 13
while true do // I. Assign vertices to cores assign up to α assignable vertices to core 1 and superstep s with Rules I & II Ω1 ← total newly assigned weight to core 1 for core p = 2, . . . , k in order do Ωp ← 0 while Ωp ̸≈ Ω1 do if ∃ assignable vertex v to core p and superstep s then assign vertex v to core p and superstep s with Rules I & II Ωp ← Ωp + ω(v) else break
18
// II. Score P assignments p Ωp β← maxp Ωp + 2000 exit ← true if the score β is high enough then deem current assignments as worthy exit ← false
19 20
if not enough assignable vertices for next superstep then exit ← true
21
undo vertex assignments with superstep s α ← 23 α if exit then redo last worthy assignments s←s+1 break inner loop
14 15 16 17
22 23 24 25 26
The algorithm proceeds by assigning vertices in the order of supersteps. During each superstep it assigns assignable vertices to cores, cf. Rule I. Rather than going through time steps to keep a work balance between cores, as you would in (barrier) list schedulers [Gra69, ACD74, HCAL89, RVG02, MSQ03, PAKY24], the algorithm has a guess α as to how long it may compute during a superstep and fills up each core up to that limit in order. This allows the algorithm to prioritise locality of the execution according to Rule II. The second nature of Rule I is that it extends compute during a single superstep by assigning more restricted vertices first. Once no further assignments are being made, the assignments in their totality
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
6
are measured by a scoring function, cf. Line 14. If the score is high enough, the assignments are considered successful and another attempt with increased α is initiated. The increase of α is exponential as to amortise the costs of retrying assignments of the current superstep. Starting with a small α and gradually increasing it allows the algorithm to keep a work balance, but also maximise independent compute. The balance between work balance and independent compute is determined by the score β. The scoring function is taken from [BPS+ 25] and is inversely proportional to the execution time of the current superstep in a bulk-synchronous-parallel model without communication [Val90a] scaled up to the whole graph. The constant 2000 thus plays the role of the synchronisation cost. Albeit our algorithm produces a stale-synchronousparallel schedule with staleness 2 rather than a bulk-synchronous-parallel schedule and can thus overlap synchronisation and compute, we find that the scoring function nevertheless does a good job in balancing parallelism and the length of independent compute. Returning to the flow of the algorithm, the superstep attempts are interrupted when either the score β starts to drop or it is deemed that there are no longer sufficient assignable vertices left for the following superstep. At this point, the best assignment measured by the score β is chosen as the final assignment of the current superstep and the algorithm moves onto the next superstep. The algorithm is efficiently implemented using sorted double-ended queues for the globally ready queues of vertices which are assignable to all cores for a given superstep, and heap for locally ready queues of vertices which are assignable to a single given core and given superstep. As a result, the presented algorithm runs in almost linear time under reasonable assumptions. We refer to [BPS+ 25, Appendix B] for a proof and §6.5 for empirical data. At last, we also remark on the choice of staleness 2. It is the smallest number that allows the overlap of compute and synchronisation. Larger numbers would only be beneficial if the compute is too small compared to the synchronisation. We find this to be unlikely the case. Moreover, in our initial testing, we have found that larger staleness negatively affects the schedule quality, thus degrading the overall performance. Hence, we chose staleness equal to 2.
4. SpTrSV implementation 4.1. Synchronisation mechanism. A stale-synchronous-parallel schedule of staleness s ≥ 2 allows one to overlap synchronisation with computation. To this end, we implemented a weak (or fuzzy) barrier [Gup89] which separates the arrive and wait operations. In our barrier implementation, each core has its own flag to signal that it has arrived at the barrier. The flag-owning core has exclusive write access to its flag. All other cores may read this flag to check whether the flag-owning core has arrived at the barrier. This allows us to avoid read-modify-write operations. In addition, our barrier makes use of caching [LBC09, Rig20], that is, every core has their own local cache of the flags of all other cores. This reduces contention and communication. To make the most out of the caching mechanism, to save memory, and to reduce code complexity, we combine all barriers into one by having a counter as the flag. This counter records the superstep which the core is currently processing or is waiting to start processing, meaning the core has completed the computation of all previous supersteps. This greatly benefits the caching mechanism as the flag of a core which is (far) ahead in terms of supersteps only needs to be checked once until the probing core has caught up. In fact, the probing core not only can catch up but also move ahead of the cached value due to the staleness before having to read again the flag of the same core. The wait operation was implemented using busy waiting. For other kinds of barriers, we refer to the survey paper [HMMR05] and the references therein.
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
7
4.2. Stale-synchronous-parallel SpTrSV kernel. Our stale-synchronous-parallel SpTrSV kernel solves a lower triangular linear system in parallel according to a stale-synchronousparallel schedule (π, σ) of staleness 2, cf. §2.2. We parallelise the kernel using OpenMP [DM98] threads. Each thread runs on its own physical core and represents a core p in the schedule. In each superstep s, every core p executes the computations associated to the rows whose corresponding vertex v has been assigned to said processor and superstep, i.e., v ∈ π −1 ({p}) ∩ σ −1 ({s}), in ascending order. The supersteps are surrounded with the weak barrier mechanics from §4.1, that is, the wait operation at the beginning of the superstep and arrive at the end of the superstep. We note that the wait operation only needs to be invoked if there is something to compute for the core during the superstep. If there is not, the wait can be safely skipped. As supersteps are sometimes empty due matrix-parallelism limitations and algorithm choices, this leads to improvements due to the implemented caching mechanism of the barrier, cf. §4.1. Algorithm 4.1: Stale-synchronous-parallel SpTrSV kernel Data: An invertible lower-triangular matrix A ∈ Rn×n , a vector b ∈ Rn , a set of cores P = {1, 2, . . . , k}, and a stale-synchronous-parallel schedule (π, σ) for GA of staleness 2, cf. §2.2. Result: A vector x ∈ Rn such that Ax = b. 1 for core p ∈ P do in parallel 2 for s = 0, 1, . . . , max σ(VA ) do 3 if π −1 ({p}) ∩ σ −1 ({s}) ̸= ∅ then 4 5
wait(p, s − 2)
for i ∈ π −1 ({p}) ∩ σ −1 ({s}) do
6 7
t ← bi for j = 1, . . . , i − 1 with Ai,j ̸= 0 do t ← t − Ai,j xj
8
xi ← t/Ai,i
9
arrive(p, s)
5. Experimental setup This section presents the benchmarks and setup used for our experiments. In large part, we follow the experimental setup of [BPSY25, BPS+ 25] and [ZCL+ 22]. The implementation of our algorithm and the test suite are available open source within the OneStopParallel scheduling framework on Github [BLM+ 24]. 5.1. Methodology. Our own scheduling algorithm is evaluated using the SpTrSV kernel described in §4.2. The kernel is parallelised through the OpenMP library [DM98], with OMP PROC BIND=close and OMP PLACES=cores as the settings. The scheduler GrowLocal is implemented in a similar fashion with the appropriate bulk-synchronous-parallel kernel implementation. The baseline schedulers HDagg and SpMP are evaluated using the Sympiler framework of Chesmi et al. [CKSD17, Che22]. This includes the original implementation of HDagg [ZCL+ 22] and has SpMP [PSSD14] already integrated into it. For every scheduler and problem instance, we repeat the SpTrSV execution 100 times in order to obtain a reliable estimation of the running time. We also add two untimed executions beforehand in order to ensure that the measurements all happen in a ‘hot’ state of the system. Between two consecutive SpTrSV executions, the right-hand-side vector b is always reset to the all-ones vector.
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
Ratio of runs within threshold
1.0
Ratio of runs within threshold
1.0
Ratio of runs within threshold
HDagg
1.0
SpMP
SuiteSparse
GrowLocal
8
ElasticDivide
AMD iCholesky
Metis ND
ARM Kunpeng
0.8 0.6 0.4 0.2 0.0
0.8
x86 EPYC
0.6 0.4 0.2 0.0
0.8
x86 Xeon
0.6 0.4 0.2 0.0 1.00
1.25
1.50
1.75
2.00
2.25
Performance threshold
2.50
2.75
3.00 1.00
1.25
1.50
1.75
2.00
2.25
Performance threshold
2.50
2.75
3.00 1.00
1.25
1.50
1.75
2.00
2.25
2.50
2.75
3.00
Performance threshold
Figure 5.1. Performance profiles of the various SpTrSV scheduling algorithms on a given data set and CPU architecture. In each plot, the x-axis represents a threshold and the y-axis the ratio of runs of a given algorithm which are within this threshold times the best run of the same matrix and core count.
The test suite of our kernel inherently implements the process above. With minimal modifications, we also adjusted the runtime measurement module of Sympiler to follow the same methodology. All scheduling algorithms (ours and the baselines) are implemented in C++ and were compiled with GCC (11.4.0 or 11.5.0) using the optimisation flag -O3. The running times of the algorithms are measured using the high-resolution clock in C++ (available in std::chrono). 5.2. Data sets. We evaluate the SpTrSV scheduling algorithms on several different data sets. The data sets all originate from the SuiteSparse Matrix Collection [DH11]. This collection contains a diverse set of matrices from a variety of applications, and it is the prominent evaluation benchmark in previous works on SpTrSV scheduling [PSSD14, ZCL+ 22, BPSY25, BPS+ 25]. Following in large part [BPSY25, BPS+ 25] and [ZCL+ 22], one of our data sets consists of matrices directly from SuiteSparse. The other two data sets consist of modified versions of the same matrices in order to align more closely with applications and previous benchmarks. An overview of the matrices in the data sets as well as some basic information, including the number of non-zero entries and the average wavefront size, may be found in §A. 5.2.1. SuiteSparse. This data set is taken from [BPSY25, BPS+ 25] and consists of the lower triangular part of real symmetric positive definite matrices from the SuiteSparse Matrix Collection [DH11]. It encompasses 33 matrices. The SpTrSV problem associated with each matrix involves at least 2 million floating point operations and has an average wavefront size of at least 44. An overview of the properties of these matrices is available in Table A.1 of the supplement.
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
9
5.2.2. SuiteSparse Eigen incomplete Cholesky (iChol). Our second data set is also taken from [BPSY25, BPS+ 25] and consists of lower triangular matrices obtained from an incomplete Cholesky decomposition on the real symmetric matrices from the SuiteSparse data set, see §5.2.1. The incomplete Cholesky decomposition was performed using the ‘IncompleteCholesky’ method in Eigen [GJ+ 10], which entailed applying its internal reordering method ‘AMDOrdering’ [ADD96]. We note that the matrix ‘bundle adj’ segmentation-faults during this process. The properties of the matrices in this data set are outlined in Table A.2 of the supplement. 5.2.3. SuiteSparse METIS (METIS). For the final data set, we use the data set considered in [ZCL+ 22]. It is once more derived from real symmetric positive definite matrices from the SuiteSparse Matrix Collection. The matrices are preprocessed using the fill-reducing method (nested dissection) of the METIS partitioner [KK98], which symmetrically permutes the original matrices. After this preprocessing step, the lower triangular part is taken. The matrices in this data set are representative of SpTrSV workloads in a Gauß–Seidel or a zero-fill-in incomplete Cholesky preconditioned conjugate gradient method for sparse symmetric solve. An overview of the matrices in this data set is available in Table A.3 of the supplement. 5.3. CPU architectures. We evaluated our SpTrSV scheduler on both x86 and ARM architectures. The concrete specifications for the three machines used in the experiments are as follows: • Huawei Kunpeng 920-4826 (Hi1620) processor (ARM, introduced 2019), with 512 GB memory, theoretical peak memory throughput of 187.7 GB/s, and two sockets with 48 cores each for a total of 96 cores; 2 NUMA domains per socket; kernel version 5.15.0; GCC version 11.4.0. • AMD EPYC 7763 processor (x86, introduced 2021), with 1024 GB memory, theoretical peak memory throughput of 204.8 GB/s, and two sockets with 64 cores each for a total of 128 cores; 1 NUMA domain per socket; kernel version 5.15.0; GCC version 11.4.0; • Intel Xeon Gold 6238T processor (x86, introduced 2019), with 192 GB memory, theoretical peak memory throughput of 140.8 GB/s, and two sockets with 22 cores each for a total of 44 cores; 1 NUMA domain per socket; kernel version 5.14.0; GCC version 11.5.0; 6. Evaluation 6.1. Runtime. To summarise our findings, we have condensed in Figure 5.1 the SpTrSV runtimes into performance profiles [DM02]. We have generated a performance profile for each data set and processor architecture, cf. §5. In each of these performance profiles, we took the best SpTrSV runtime for each core count and matrix combination. For a given algorithm, we then computed the ratio of all SpTrSV from that algorithm that are with a given threshold of the respective best SpTrSV runtime. The resulting curve is then plotted with the x-axis representing the threshold and the y-axis representing the ratio. The closer a curve of a given algorithm is to the top left corner, the better the algorithm is at performing parallel SpTrSV. We can clearly see that the new algorithm ElasticDivide performs exceptionally well on the ARM architecture. Similarly, on x86, ElasticDivide is the most performant, though it shares the top spot with GrowLocal. We explain this with the increased NUMA domains the ARM Kunpeng architecture has, favouring asynchronous compute. SpMP is a strong contender on the SuiteSparse data set and performs well overall, but ultimately falls a bit short of ElasticDivide and GrowLocal. HDagg struggles with the hard to parallelise SuiteSparse data set and only seems to do well on x86 Xeon architecture with smaller core count.
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE HDagg
SpMP
GrowLocal
ElasticDivide
AMD iCholesky
Metis ND
30
ARM Kunpeng
Speed-up over Serial
SuiteSparse
10
20 10 0
16
32
48
64
80
96
16
32
48
64
80
96
16
32
48
64
80
96
15 10 5 0
16
32
48
64
80
96
112
128 16
32
48
64
80
96
112
128 16
32
48
64
80
96
112
128
25 20
x86 Xeon
Speed-up over Serial
x86 EPYC
Speed-up over Serial
20
15 10 5 0
11
22
33
Cores
44
11
22
33
Cores
44
11
22
33
44
Cores
Figure 6.1. Scaling plots depicting the geometric mean speed-up over Serial and interquartile range of each algorithm on a given data set, architecture, and core count. 6.2. Scaling. We investigate how the algorithms scale with the core count by graphing the geometric mean and interquartile range of the speed-up over serial execution. We do this in Figure 6.1 for every pair of architecture and data set, cf. §5. We once more see that ElasticDivide outperforms all other algorithms on all data sets on the ARM architecture. It reaches a geometric mean speed-up of 20.00 on the Metis data set on 64 cores. This is 12% higher than the speed-up of 17.93 reached by the second best algorithm GrowLocal on the same data set (and architecture). On the x86 architecture, GrowLocal and ElasticDivide compete closely for the best algorithm. However, as the core count increases, in particular on the SuiteSparse data set, SpMP continues to improve relative to GrowLocal and ElasticDivide and in some cases matches or improves upon their performance. HDagg does not scale well with the number of cores except for the Metis data set as the nested dissecting reordering works well with the wavefront-gluing method that HDagg employs. On all architectures, we see a drop in speed-ups at the 80 core count. On the EPYC processor this coincides with a second socket, which is a likely contributor, but on the Kunpeng processor we already require the second socket at 64 cores. Thus, we find the more likely explanation is that the parallelism which is exposed by the matrices and then picked up by the algorithms peters out around that core count, at least for a majority of matrices. The reader may consult the Tables A.1, A.2, and A.3 for the average wavefront size of the matrices. 6.3. Synchronisation. In Table 6.1, we present the geometric-mean reduction of the number of synchronisation barriers (weak or strong) relative to the number of wavefronts of the matrix. In other words, the reduction of synchronisation barriers compared to a wavefront scheduler. The numbers we present are from all data sets over a sample of core counts. We note that SpMP is missing in the table as it employs an asynchronous (point-to-point) method of compute.
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
Data
Suite
iChol
Metis
Cores
HDagg
GrowLocal
ElasticDivide
16
1.34
16.02
6.74
32
1.15
12.78
5.10
64
1.06
11.35
4.48
128
1.07
10.69
4.17
16
1.64
20.74
10.33
32
1.52
17.35
8.62
64
1.44
14.58
7.76
128
1.37
13.30
7.01
16
2.57
16.73
8.59
32
2.10
14.54
8.15
64
1.87
12.70
7.22
128
1.74
11.38
6.73
11
Table 6.1. Geometric-mean reduction of the number of synchronisation barriers over the number of wavefronts for each algorithm. We see that the wavefront-gluing technique of HDagg is more effective on the Metis nested dissection reordered data set and also at lower core counts. At higher core counts, in particular on the SuiteSparse data set, it struggles to glue together wavefronts. Under these settings, HDagg only becomes marginally better than a wavefront scheduler, which perhaps explains its lack of performance in the higher core count numbers. Compared to the breadth-first search of wavefront-based schedulers like HDagg, the depthfirst search approaches of GrowLocal and ElasticDivide are able to pack in more compute before having to synchronise. Albeit the reported reduction by ElasticDivide is about half the one by GrowLocal, one has to recall that the synchronisations by ElasticDivide are weak and since the staleness is 2, cf. §2.2, its numbers should be multiplied by two. Having taken that into account, we see that ElasticDivide does slightly better than GrowLocal on the Metis data set and slightly worse on the SuiteSparse data set. 6.4. Overlap of Synchronisation and Compute. In order to investigate the effect of overlapping synchronisation and compute, we execute the stale-synchronous-parallel schedule computed by ElasticDivide twice: once using our stale-synchronous-parallel SpTrSV kernel, cf. §4.2, and once using a bulk-synchronous-parallel SpTrSV kernel, implemented using OpenMP [DM98]. We note that this results in a valid execution, cf. Remark 2.3. In Table 6.2, we display the relative speed-up that the stale-synchronous-parallel SpTrSV kernel achieves over the bulk-synchronous-parallel one. The data shows a clear picture: overlapping synchronisation and compute is beneficial. In general, we see higher improvements for data sets whose matrices are harder to parallelise such as the SuiteSparse data set and lower numbers for easier to parallelise matrices which are present in the Metis nested dissection data set. We also see that the costlier the synchronisation is, for example through higher non-uniform memory architectures, the larger the improvement is. That is, we see larger speed-ups in ARM than x86 and larger speed-ups for higher core count. There is one exception to the latter and that is the x86 EPYC architecture on the Metis data set. This phenomenon requires further investigation, though we offer possible explanations: the compiler produces a more efficient type of barrier for this architecture than our barrier, cf. §4.1 and [HMMR05], and/or the stale-synchronous-parallel schedule looks kind of like a bulk-synchronous-parallel schedule, that is every other superstep is almost empty.
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
Architecture
ARM Kunpeng
x86 EPYC
x86 Xeon
Cores
Suite
iChol
Metis
16
1.57
1.07
1.06
32
1.95
1.15
1.13
48
2.17
1.35
1.21
64
2.82
1.37
1.34
80
3.23
1.76
1.48
96
4.51
2.02
1.65
16
1.23
1.08
1.07
32
1.36
1.09
1.07
48
1.47
1.12
1.06
64
1.60
1.12
1.05
80
1.84
1.13
1.03
96
1.93
1.18
1.01
112
2.01
1.16
0.96
128
2.19
1.13
0.99
11
1.11
1.05
1.03
22
1.33
1.08
1.12
33
1.34
1.13
1.12
44
1.15
1.10
1.08
12
Table 6.2. Geometric-mean speed-ups of stale-synchronous-parallel execution over bulk-synchronous-parallel execution for all architectures, core counts, and data sets.
HDagg
SpMP
ARM Kunpeng
GrowLocal
ElasticDivide
x86 EPYC
x86 Xeon
Schedule compute time [s]
103 102 101 100 10−1 10−2 106
107
Number of non-zeroes
108
106
107
Number of non-zeroes
108
106
107
108
Number of non-zeroes
Figure 6.2. Scheduling times of each algorithm plotted against number of non-zeroes of a matrix. The best ℓ2 -fitted line of the shape log(y) = m·log(x)+c is also depicted.
6.5. Amortisation. An important aspect of the parallel SpTrSV schedules is the time it takes to generates them in the first place and how the quality of the schedule compares to the overhead to generate it in the first place. To this end, we look at how the scheduling time increases with the number of non-zeroes in the matrix. We display the results in Figure 6.2. We also compute the amortisation costs associated with a schedule [ZCL+ 22, BPSY25, BPS+ 25]. This is defined as the number of executions at which point it is more beneficial to compute a parallel schedule and then execute in parallel, opposed to always using serial execution:
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
Architecture
Data
ARM Kunpeng
x86 EPYC
x86 Xeon
13
Algorithm HDagg
SpMP
GrowLocal
ElasticDivide
Metis
54.82
—
33.59
13.60
iChol
211.01
—
35.38
13.16
Suite
1763.03
—
31.10
11.59
Metis
85.59
5.41
52.91
21.32
iChol
331.20
7.13
53.88
22.97
Suite
2819.82
9.20
58.28
17.51
Metis
59.52
5.88
43.79
21.60
iChol
227.06
7.54
46.66
22.80
Suite
1398.14
7.20
45.51
25.10
Table 6.3. Median amortisation costs of each algorithm on a full single socket for each architecture and data set.
( namortisation =
tschedule , tserial −tparallel
+∞,
if tserial > tparallel , otherwise,
(6.1)
where the bar on top of the times denotes arithmetic mean. In Table 6.3, we record the median amortisation costs for each algorithm, data set, and architecture, cf. §5. As the core count, we have taken the number of cores corresponding to a single socket. It is apparent that SpMP is the fastest in scheduling time and that it scales linearly with the number of non-zeroes. This goes to show the engineering efforts that went into the algorithm, in particular the efforts to the make the scheduling algorithm parallel in the first place. ElasticDivide and GrowLocal generate their schedule using a serial algorithm and are therefore slower than SpMP. Nevertheless, their runtime scales (provably) linearly with the number of non-zeroes. HDagg is considerably slower and does not scale linearly. Indeed, the slopes in Figure 6.2 of the best ℓ2 -fitted lines of the shape log(y) = m · log(x) + c are between 1.38 and 1.56 for HDagg, whereas they are ≤ 1.0 for all other algorithms. Albeit that ElasticDivide and GrowLocal follow a similar algorithm design, the distance between their ℓ2 -fitted lines in Figure 6.2 indicates that ElasticDivide is twice as fast in computing a schedule. We attribute this to a better implementation with improved data structures. The amortisation cost table, Table 6.3, reflects the scheduling time data depicted in Figure 6.2. SpMP comes in the lowest with an amortisation cost of 5-10 across x86 architectures. ElasticDivide achieves low amortisation costs, ranging from 12 to 25 across all architectures and data sets. GrowLocal has costs ranging from 31 to 58, performing best on the Metis data set and worst on the SuiteSparse data set. HDagg has the largest amortisation costs, ranging from 55-86 on the most favourable Metis data set all the way up to 2820 on the SuiteSparse data set on x86 EPYC. 7. Conclusion Synchronising many cores is expensive, more so in non-uniform memory architectures. To combat this expense, one may reduce the number of synchronisations necessary through preprocessing the matrix whenever possible, e.g., nested dissecting, through algorithms that reduce the number of synchronisations to a minimum, and, as we have shown in this paper, through overlapping compute with synchronisation with execution models such as stalesynchronous-parallel. The latter can improve performance easily by 10-30% and when the
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
14
number of synchronisations is large and the architecture is highly non-uniform by 100% and in some cases up to 350%, cf. §6.4. Our stale-synchronous-parallel scheduler, ElasticDivide, demonstrates that the fine-grained dependency graphs present in sparse triangular linear systems allow for the necessary flexibility to fit into a stale-synchronous-parallel schedule without (significant) increase in the number of synchronisation events when appropriately normalised. This resulted in geometric-mean speed-ups of around 7-30% against state-of-the-art scheduler GrowLocal on ARM Kunpeng, with even higher speed-ups for higher core and NUMA-domain count, cf. §6.2. Albeit our presentation and evaluation was conducted on a variety of CPUs, our algorithm ElasticDivide and more generally our ideas and findings are transferable to other compute architectures. For instance, ElasticDivide is directly applicable to solving sparse triangular systems on accelerated compute systems such as GPUs. There it could be used to mask CPU to GPU and cross-GPU latencies with compute. Appendix A. Tables of matrices The tables provided here give more details on the the matrices used in the experiments, cf. §5.2, together with some basic statistics. Matrix
Size
#Non-zeroes
Average wavefront size
af 0 k101
503,625
9,027,150
74
af shell7
504,855
9,046,865
135
apache2
715,176
2,766,523
1,077
audikw 1
943,695
39,297,771
203
bmw7st 1
141,347
3,740,507
199
bmwcra 1
148,770
5,396,386
204
bone010
986,703
36,326,514
470
boneS01
127,224
3,421,188
156
boneS10
914,898
28,191,660
386
Bump 2911
2,911,419
65,320,659
283
bundle adj
513,351
10,360,701
57,039
consph
83,334
3,046,907
139
Dubcova3
146,689
1,891,669
44
ecology2
999,999
2,997,995
500
Emilia 923
923,136
20,964,171
176
Fault 639
638,802
14,626,683
143
Flan 1565
1,564,794
59,485,419
200
G3 circuit
1,585,478
4,623,152
611
Geo 1438
1,437,960
32,297,325
246
220,542
5,494,489
365
1,498,023
31,207,734
95
503,712
18,660,027
287
hood Hook 1498 inline 1
Table A.1. Matrices used in the evaluation from the SuiteSparse Matrix Collection [DH11]. The average wavefront size has been rounded down. (Continued on next page)
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
15
Matrix
Size
#Non-zeroes
Average wavefront size
ldoor
952,203
23,737,339
141
msdoor
415,863
10,328,399
59
offshore
259,789
2,251,231
75
parabolic fem
525,825
2,100,225
75,117
PFlow 742
742,793
18,940,627
118
Queen 4147
4,147,110
166,823,197
342
s3dkt3m2
90,449
1,921,955
60
Serena
1,391,349
32,961,525
298
shipsec1
140,874
3,977,139
67
StocF-1465
1,465,137
11,235,263
487
thermal2
1,228,045
4,904,179
991
Table A.1. Matrices used in the evaluation from the SuiteSparse Matrix Collection [DH11]. The average wavefront size has been rounded down.
Matrix
Size
#Non-zeroes
Average wavefront size
af 0 k101 iCh
503,625
9,027,150
195
af shell7 iCh
504,855
9,046,865
668
apache2 iCh
715,176
2,766,523
79,464
audikw 1 iCh
943,695
39,297,771
138
bmw7st 1 iCh
141,347
3,740,507
340
bmwcra 1 iCh
148,770
5,396,386
89
bone010 iCh
986,703
36,326,514
340
boneS01 iCh
127,224
3,421,188
245
boneS10 iCh
914,898
28,191,660
521
Bump 2911 iCh
2,911,419
65,320,659
1,048
consph iCh
83,334
3,046,907
78
Dubcova3 iCh
146,689
1,891,669
1,594
ecology2 iCh
999,999
2,997,995
142,857
Emilia 923 iCh
923,136
20,964,171
511
Fault 639 iCh
638,802
14,626,683
422
Flan 1565 iCh
1,564,794
59,485,419
689
G3 circuit iCh
1,585,478
4,623,152
88,082
Geo 1438 iCh
1,437,960
32,297,325
768
220,542
5,494,489
1,050
hood iCh
Table A.2. Matrices used in the evaluation from SuiteSparse Matrix Collection [DH11]. These matrices were transformed using the incomplete Cholesky method with AMD reordering of Eigen [GJ+ 10]. The average wavefront size has been rounded down. (Continued on next page)
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
Matrix
16
Size
#Non-zeroes
Average wavefront size
1,498,023
31,207,734
649
inline 1 iCh
503,712
18,660,027
679
ldoor iCh
952,203
23,737,339
3,317
msdoor iCh
415,863
10,328,399
956
offshore iCh
259,789
2,251,231
1,114
parabolic fem iCh
525,825
2,100,225
19,475
PFlow 742 iCh
742,793
18,940,627
240
Queen 4147 iCh
4,147,110
166,823,197
719
s3dkt3m2 iCh
90,449
1,921,955
104
Serena iCh
1,391,349
32,961,525
940
shipsec1 iCh
140,874
3,977,139
259
StocF-1465 iCh
1,465,137
11,235,263
2,990
thermal2 iCh
1,228,045
4,904,179
47,232
Hook 1498 iCh
Table A.2. Matrices used in the evaluation from SuiteSparse Matrix Collection [DH11]. These matrices were transformed using the incomplete Cholesky method with AMD reordering of Eigen [GJ+ 10]. The average wavefront size has been rounded down.
Matrix
Size
#Non-zeroes
Average wavefront size
af 0 k101 metis
503,625
9,027,150
610
af shell10 metis
1,508,065
27,090,195
1,065
apache2 metis
715,176
2,766,523
47,678
audikw 1 metis
943,695
39,297,771
1,734
bmwcra 1 metis
148,770
5,396,386
473
bone010 metis
986,703
36,326,514
1,326
boneS10 metis
914,898
28,191,660
2,401
bundle adj metis
513,351
10,360,701
11,407
cant metis
62,451
2,034,917
333
consph metis
83,334
3,046,907
247
crankseg 2 metis
63,838
7,106,348
86
ecology2 metis
999,999
2,997,995
62,499
Emilia 923 metis
923,136
20,964,171
2,107
Fault 639 metis
638,802
14,626,683
1,458
Flan 1565 metis
1,564,794
59,485,419
2,569
Table A.3. Matrices used in the evaluation from the SuiteSparse Matrix Collection [DH11]. These matrices were symmetrically permuted using the fill-reducing method ‘METIS NodeND’ of [KK98]. The average wavefront size has been rounded down. (Continued on next page)
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
Matrix
17
Size
#Non-zeroes
Average wavefront size
G3 circuit metis
1,585,478
4,623,152
93,263
Geo 1438 metis
1,437,960
32,297,325
2,887
gyro metis
17,361
519,260
88
hood metis
220,542
5,494,489
984
1,498,023
31,207,734
4,059
inline 1 metis
503,712
18,660,027
1,549
ldoor metis
952,203
23,737,339
4,858
m t1 metis
97,578
4,925,574
268
msdoor metis
415,863
10,328,399
1,856
nasasrb metis
54,870
1,366,097
287
PFlow 742 metis
742,793
18,940,627
1,023
pwtk metis
217,918
5,926,171
511
raefsky4 metis
19,779
674,195
111
ship 003 metis
121,728
4,103,881
494
shipsec8 metis
114,919
3,384,159
456
StocF-1465 metis
1,465,137
11,235,263
11,446
thermal2 metis
1,228,045
4,904,179
45,483
tmt sym metis
726,713
2,903,837
26,915
x104 metis
108,384
5,138,004
306
Hook 1498 metis
Table A.3. Matrices used in the evaluation from the SuiteSparse Matrix Collection [DH11]. These matrices were symmetrically permuted using the fill-reducing method ‘METIS NodeND’ of [KK98]. The average wavefront size has been rounded down. Appendix B. Flops/s In Figures B.1, B.2, and B.3, we depict the double precision floating point operations per second on each graph, data set, and architecture. The number of cores was chosen as a single socket. Error bars indicate the standard deviation of the measurements taken.
GFP64/s
0.0
2.5
5.0
7.5
10.0
12.5
15.0
17.5
20.0
0
20
40
60
80
0
5
10
15
20
100
GFP64/s
she
ll7
che
2
ikw
w7 st wc ra
1
1
1
ne0
10
neS neS
mp
01 10
291 1
nd
le
adj sph
cov
a3 HDagg
y2
ilia
923 SpMP
639
n1
565 GrowLocal
circ
uit
o1
438 od
ElasticDivide
ok
149
8
or
ldo
1
ne
ore
abo
lic
fem
low
een
742
414
7
kt3
m2
ena
pse
c1
cF-
146 5 al2 x86 Xeon
Figure B.1. Double precision floating point operations per second on the SuiteSparse data set using number of cores equal to one socket.
GFP64/s
25
0k 101
30
af
af
apa aud bm
bm
bo bo bo Bu
bu
con Du b
eco log Em
Fau lt Fla G3 Ge
ho Ho
inli
or
do
ms
off sh par
PF Qu
s3d
Ser shi Sto
rm
x86 EPYC
the
ARM Kunpeng
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
18
0.0
2.5
5.0
7.5
10.0
12.5
15.0
0
10
20
30
40
50
60
70
0
5
10
17.5
GFP64/s 15
20.0
GFP64/s
Ch
po
ol
stC
po
hol
stC h
ol
ost
Ch
stC
ol
hol ol
stC
hol
stC h
ol
ost
Ch
ol
ost
Ch
ol
ost
Ch
ol
ost
Ch
ol
stC
HDagg
hol
po
stC
hol
po
stC h
po
stC
SpMP
ol
hol
stC
GrowLocal
od
hol
ost
Ch
ol
po
stC hol ElasticDivide
stC
hol
stC
hol
ost
Ch
ol
ost
Ch
ol
po
stC h
po
ol
stC
742
hol
po
stC
ee
7p
hol
ost
Ch
k
ol
po
stC
hol
po
stC
hol
po
stC hol
cF
po
stC
hol
po
stC h
ol x86 Xeon
Figure B.2. Double precision floating point operations per second on the AMD reordered incomplete Cholesky data set using number of cores equal to one socket.
GFP64/s
20
ost
25
1p
e2
1p
w
po
wc r
ost Ch
po po
0p
1p
hp 3p
po
ili
39
n
156 5
c
po
8p
ok
po
po
or p or p
sho re
li
cf em
low
30
af
0k 10
af
she ll7
apa ch au
dik w
bm
7st 1
bm
a1 p
b
one 010
bo
neS 01
b
one S1
Bu m
p2 91 con sp
Du
bco va
eco
log y2
Em
a9 23
Fa
ult 6
Fla G3
ircu it
Ge o1 43 ho Ho
149 8
in
line 1 ldo
m
sdo
off par ab o PF Qu
n4 14
s3d
t3m 2
Ser ena shi
pse c1
Sto
-14 65
th
l2
x86 EPYC
erm a
ARM Kunpeng
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
19
0.0
2.5
5.0
7.5
10.0
12.5
15.0
0
10
20
30
40
50
60
70
0
5
10
17.5
GFP64/s 15
20.0
GFP64/s
me
tis
me ti
s
che
2m eti
s
ikw
1m
eti
s
1m eti
s
ne0
10
10
me ti
s
me ti
s
dle
adj
me ti
s
tm
sph
eti
s
me
tis
kse
g2
me ti
s
log
y2 me ti
s
923
me ti
s
t6
me
HDagg
39
tis
565
me
tis
ircu
SpMP
it m eti
s
o1
438
me
tis
om eti
me ti
GrowLocal
s
od
s
149
8m eti
s
1m
ElasticDivide
ne
eti
s
or m
eti
t1
s
me
tis
or m
eti
asr b
s
me
tis
742
me ti
s
me
tis
me ti
s
p0
03
me
tis
pse
c8
me ti
s
146
5m
eti
rm al2
s
me
tis
sym
me
tis
4m
eti
s x86 EPYC
x86 Xeon
Figure B.3. Double precision floating point operations per second on the Metis data set using number of cores equal to one socket.
GFP64/s
20
101
25
0k
ll10
30
af a
fs he apa
aud b
mw cra bo b
one S
bu n
can con cra n
eco
Em ilia Fau l
Fla n1 G3 c Ge
gyr ho Ho ok
inli
ldo m
ms do nas PF low
pw tk rae fsk y4 shi
shi S
toc Fthe
tm t
x10
ARM Kunpeng
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
20
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
21
References [ACD74]
Thomas L. Adam, K. Mani Chandy, and J. R. Dickson. A comparison of list schedules for parallel processing systems. Communications of the ACM, 17(12):685–690, 1974. [ADD96] Patrick R. Amestoy, Timothy A. Davis, and Iain S. Duff. An approximate minimum degree ordering algorithm. SIAM Journal on Matrix Analysis and Applications, 17(4):886–905, 1996. [APc04] Cevdet Aykanat, Ali Pinar, and Ümit V. Çatalyürek. Permuting sparse rectangular matrices into block-diagonal form. SIAM Journal on Scientific Computing, 25(6):1860–1879, 2004. [AS89] Edward Anderson and Youcef Saad. Solving sparse triangular linear systems on parallel computers. International Journal of High Speed Computing, 1(01):73–95, 1989. [AYU21] Najeeb Ahmad, Buse Yilmaz, and Didem Unat. A split execution model for sptrsv. IEEE Transactions on Parallel and Distributed Systems, 32(11):2809–2822, 2021. [BAvL+ 19] Rob H. Bisseling, Bas Fagginger Auer, Tristan van Leeuwen, Wouter Meesen, Marco van Oort, Daan Pelt, Brendan Vastenhouw, and Albert-Jan N. Yzelman. Mondriaan (version 4.2.1). http: //www.staff.science.uu.nl/~bisse101/Mondriaan/mondriaan_v4.2.1.tar.gz, 2019. [BLM+ 24] Toni Böhnlein, Benjamin Lozes, Christos K. Matzoros, Pál András Papp, and Raphael S. Steiner. OneStopParallel. https://github.com/Algebraic-Programming/OneStopParallel, 2024. [BPS+ 25] Toni Böhnlein, Pál András Papp, Raphael S. Steiner, Christos K. Matzoros, and AlbertJan N. Yzelman. Efficient parallel scheduling for sparse triangular solvers. arXiv preprint arXiv:2503.05408, 2025. [BPSY25] Toni Böhnlein, Pál András Papp, Raphael S. Steiner, and Albert-Jan N. Yzelman. Efficient parallel scheduling for sparse triangular solvers. In 2025 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW), pages 1263–1265, 2025. [CAGL+ 87] C. Cleveland Ashcraft, Roger G. Grimes, John G. Lewis, Barry W. Peyton, Horst D. Simon, and Petter E. Bjørstad. Progress in sparse matrix methods for large linear systems on vector supercomputers. The International Journal of Supercomputing Applications, 1(4):10–30, 1987. [CCH+ 14] Henggang Cui, James Cipar, Qirong Ho, Jin Kyu Kim, Seunghak Lee, Abhimanu Kumar, Jinliang Wei, Wei Dai, Gregory R. Ganger, Phillip B. Gibbons, Garth A. Gibson, and Eric P. Xing. Exploiting bounded staleness to speed up big data analytics. In 2014 USENIX Annual Technical Conference (USENIX ATC 14), pages 37–48, Philadelphia, PA, June 2014. USENIX Association. [Che22] Kazem Cheshmi. Transforming Sparse Matrix Computations. PhD thesis, University of Toronto, Computer Science, 2022. [CHK+ 13] James Cipar, Qirong Ho, Jin Kyu Kim, Seunghak Lee, Gregory R. Ganger, Garth Gibson, Kimberly Keeton, and Eric P. Xing. Solving the straggler problem with bounded staleness. In 14th Workshop on Hot Topics in Operating Systems (HotOS XIV), Santa Ana Pueblo, NM, May 2013. USENIX Association. [CKSD17] Kazem Cheshmi, Shoaib Kamil, Michelle Mills Strout, and Maryam Mehri Dehnavi. Sympiler: Transforming sparse matrix codes by decoupling symbolic analysis. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC ’17, pages 13:1–13:13, New York, NY, USA, 2017. ACM. [DH11] Timothy A. Davis and Yifan Hu. The University of Florida sparse matrix collection. ACM Transactions on Mathematical Software (TOMS), 38(1):1–25, 2011. [DM98] Leonardo Dagum and Ramesh Menon. Openmp: an industry standard api for shared-memory programming. Computational Science & Engineering, IEEE, 5(1):46–55, 1998. [DM02] Elizabeth D. Dolan and Jorge J. Moré. Benchmarking optimization software with performance profiles. Mathematical programming, 91:201–213, 2002. [DPY25] Kiril Dichev, Filip Pawlowski, and Albert-Jan N. Yzelman. Faster distributed inference-only recommender systems via bounded lag synchronous collectives. arXiv preprint arXiv:2512.19342, 2025. [DS05] Iain S. Duff and Jennifer A. Scott. Stabilized bordered block diagonal forms for parallel sparse solvers. Parallel Computing, 31(3):275–289, 2005. [FZW+ 23] Xu Fu, Bingbin Zhang, Tengcheng Wang, Wenhao Li, Yuechen Lu, Enxin Yi, Jianqi Zhao, Xiaohan Geng, Fangying Li, Jingwen Zhang, Zhou Jin, and Weifeng Liu. Pangulu: A scalable regular two-dimensional block-cyclic sparse direct solver on distributed heterogeneous systems. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC ’23, New York, NY, USA, 2023. Association for Computing Machinery. [GBDD10] Laura Grigori, Erik G. Boman, Simplice Donfack, and Timothy A. Davis. Hypergraph-based unsymmetric nested dissection ordering for sparse lu factorization. SIAM Journal on Scientific Computing, 32(6):3426–3446, 2010. [Geo73] Alan George. Nested dissection of a regular finite element mesh. SIAM Journal on Numerical Analysis, 10(2):345–363, 1973. [GJ+ 10] Gaël Guennebaud, Benoı̂t Jacob, et al. Eigen v3. http://eigen.tuxfamily.org, 2010.
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
[Gra69]
22
Ronald L. Graham. Bounds on multiprocessing timing anomalies. SIAM journal on Applied Mathematics, 17(2):416–429, 1969. [Gup89] Rajiv Gupta. The fuzzy barrier: a mechanism for high speed synchronization of processors. In Proceedings of the Third International Conference on Architectural Support for Programming Languages and Operating Systems, ASPLOS III, page 54–63, New York, NY, USA, 1989. Association for Computing Machinery. [HCAL89] Jing-Jang Hwang, Yuan-Chieh Chow, Frank D. Anger, and Chung-Yee Lee. Scheduling precedence graphs in systems with interprocessor communication times. siam journal on computing, 18(2):244–257, 1989. [HCC+ 13] Qirong Ho, James Cipar, Henggang Cui, Jin Kyu Kim, Seunghak Lee, Phillip B. Gibbons, Garth A. Gibson, Gregory R. Ganger, and Eric P. Xing. More effective distributed ml via a stale synchronous parallel parameter server. In Proceedings of the 27th International Conference on Neural Information Processing Systems - Volume 1, NIPS’13, page 1223–1231, Red Hook, NY, USA, 2013. Curran Associates Inc. [HMMR05] Torsten Hoefler, Torsten Mehlan, Frank Mietke, and Wolfgang Rehm. A survey of barrier algorithms for coarse grained supercomputers. Universitätsbibliothek Chemnitz, 2005. [KK98] George Karypis and Vipin Kumar. A fast and high quality multilevel scheme for partitioning irregular graphs. SIAM Journal on scientific Computing, 20(1):359–392, 1998. [LBC09] Patrick P. C. Lee, Tian Bu, and Girish Chandranmenon. A lock-free, cache-efficient shared ring buffer for multi-core architectures. In Proceedings of the 5th ACM/IEEE Symposium on Architectures for Networking and Communications Systems, ANCS ’09, page 78–79, New York, NY, USA, 2009. Association for Computing Machinery. [LDS+ 23] Yang Liu, Nan Ding, Piyush Sao, Samuel Williams, and Xiaoye Sherry Li. Unified communication optimization strategies for sparse triangular solver on cpu and gpu clusters. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC ’23, New York, NY, USA, 2023. Association for Computing Machinery. [Li05] Xiaoye S. Li. An overview of superlu: Algorithms, implementation, and user interface. ACM Trans. Math. Softw., 31(3):302–325, September 2005. [LLH+ 16] Weifeng Liu, Ang Li, Jonathan Hogg, Iain S. Duff, and Brian Vinter. A synchronization-free algorithm for parallel sparse triangular solves. In Euro-Par 2016: Parallel Processing: 22nd International Conference on Parallel and Distributed Computing, Grenoble, France, August 24-26, 2016, Proceedings 22, pages 617–630. Springer, 2016. [LNL20] Zhengyang Lu, Yuyao Niu, and Weifeng Liu. Efficient block algorithms for parallel sparse triangular solve. In Proceedings of the 49th International Conference on Parallel Processing, pages 1–11, 2020. [LNP93] Joseph W. H. Liu, Esmond G. Ng, and Barry W. Peyton. On finding supernodes for sparse matrix computations. SIAM Journal on Matrix Analysis and Applications, 14(1):242–252, 1993. [LRT79] Richard J. Lipton, Donald J. Rose, and Robert Endre Tarjan. Generalized nested dissection. SIAM Journal on Numerical Analysis, 16(2):346–358, 1979. [May09] Jan Mayer. Parallel algorithms for solving linear systems with sparse triangular matrices. Computing, 86:291–312, 2009. [MSQ03] Shang Mingsheng, Sun Shixin, and Wang Qingxian. An efficient parallel scheduling algorithm of dependent task graphs. In Proceedings of the Fourth International Conference on Parallel and Distributed Computing, Applications and Technologies, pages 595–598. IEEE, 2003. [PAKY24] Pál András Papp, Georg Anegg, Aikaterini Karanasiou, and Albert-Jan N. Yzelman. Efficient Multi-Processor Scheduling in Increasingly Realistic Models. In Proceedings of the 36th ACM Symposium on Parallelism in Algorithms and Architectures. ACM, 2024. [PSSD14] Jongsoo Park, Mikhail Smelyanskiy, Narayanan Sundaram, and Pradeep Dubey. Sparsifying synchronization for high-performance shared-memory sparse triangular solver. In Supercomputing: 29th International Conference, ISC 2014, Leipzig, Germany, June 22-26, 2014. Proceedings 29, pages 124–140. Springer, 2014. [RG92] Edward Rothberg and Anoop Gupta. Parallel ICCG on a hierarchical memory multiprocessor—addressing the triangular solve bottleneck. Parallel Computing, 18(7):719–741, 1992. [Rig20] Erik Rigtorp. SPSCQueue. https://github.com/rigtorp/SPSCQueue, 2020. [RVG02] Andrei Radulescu and Arjan J. C. Van Gemund. Low-cost task scheduling for distributed-memory machines. IEEE Transactions on Parallel and Distributed Systems, 13(6):648–658, 2002. [Sal90] Joel H. Saltz. Aggregation methods for solving sparse triangular systems on multiprocessors. SIAM journal on scientific and statistical computing, 11(1):123–144, 1990. [SG04] Olaf Schenk and Klaus Gärtner. Solving unsymmetric sparse systems of linear equations with pardiso. Future Generation Computer Systems, 20(3):475–487, 2004. Selected numerical algorithms. [SGFS01] Olaf Schenk, Klaus Gärtner, Wolfgang Fichtner, and Andreas Stricker. Pardiso: a highperformance serial and parallel sparse linear solver in semiconductor device simulation. Future
ELASTICITY IN PARALLEL SPARSE TRIANGULAR SOLVE
[SMB88]
[Val90a] [Val90b] [YRE20]
[YSAU20]
[ZALC19]
[ZCL+ 22]
23
Generation Computer Systems, 18(1):69–78, 2001. I. High Performance Numerical Methods and Applications. II. Performance Data Mining: Automated Diagnosis, Adaption, and Optimization. Joel H. Saltz, Ravi Mirchandaney, and Doug Baxter. Run-time parallelization and scheduling of loops. Technical report, Institute for Computer Applications in Science and Engineering, NASA Langley Research Center, 1988. Leslie G. Valiant. A bridging model for parallel computation. Communications of the ACM, 33(8):103–111, 1990. Leslie G. Valiant. General purpose parallel architectures. In Algorithms and Complexity, pages 943–971. Elsevier, 1990. Ichitaro Yamazaki, Sivasankaran Rajamanickam, and Nathan Ellingwood. Performance portable supernode-based sparse triangular solver for manycore architectures. In Proceedings of the 49th International Conference on Parallel Processing, pages 1–11, 2020. Buse Yılmaz, Buğrra Sipahioğrlu, Najeeb Ahmad, and Didem Unat. Adaptive level binning: A new algorithm for solving sparse triangular systems. In Proceedings of the International Conference on High Performance Computing in Asia-Pacific Region, pages 188–198, 2020. Xing Zhao, Aijun An, Junfeng Liu, and Bao Xin Chen. Dynamic stale synchronous parallel distributed training for deep learning. In 39th IEEE International Conference on Distributed Computing Systems, ICDCS 2019, Dallas, TX, USA, July 7-10, 2019, pages 1507–1517. IEEE, 2019. Behrooz Zarebavani, Kazem Cheshmi, Bangtian Liu, Michelle Mills Strout, and Maryam Mehri Dehnavi. HDagg: hybrid aggregation of loop-carried dependence iterations in sparse matrix computations. In 2022 IEEE International Parallel and Distributed Processing Symposium (IPDPS), pages 1217–1227. IEEE, 2022.
Huawei, Thurgauerstrasse 80, 8050 Zurich, CH Email address: [email protected] Huawei, Thurgauerstrasse 80, 8050 Zurich, CH Email address: [email protected] Huawei, Thurgauerstrasse 80, 8050 Zurich, CH Email address: [email protected] Huawei, Thurgauerstrasse 80, 8050 Zurich, CH Email address: [email protected] Huawei, Thurgauerstrasse 80, 8050 Zurich, CH Email address: [email protected]