An Exact Parallel Dense Gauss-Seidel Algorithm with Trade-Offs between Parallel Execution Time and Preprocessing Time

An Exact Parallel Dense Gauss-Seidel Algorithm with Trade-Offs between Parallel Execution Time and Preprocessing Time⋆

arXiv:2610.39990v1 [cs.DC] 29 Oct 2026

Mugurel-Ionuţ Andreicaa , Pantelimon George Popescua , Eliana-Dina Tîrşaa a Computer Science Department Politehnica University of Bucharest Splaiul Independenţei 313, sector 6, 060042, Bucharest, Romania

Abstract The Gauss-Seidel iterative method is commonly used for solving systems of linear equations because of its good convergence properties. However, the Gauss-Seidel algorithm cannot be easily parallelized because of the data dependencies existing at each iteration. The best parallel execution time for dense systems of equations is O(n) given that enough processors are available and no preprocessing is performed. In this paper we propose a new parallel Gauss-Seidel algorithm which is capable of achieving better parallel runtimes at the expense of more preprocessing time. The algorithm can be tuned in order to provide convenient trade-offs between the preprocessing time and the best parallel runtime achievable by the algorithm. 1. Introduction In order to solve large systems of linear equations, iterative methods like Jacobi or Gauss-Seidel [1] are used. These methods produce good approximate results much faster than solving the system of equations exactly. When implemented sequentially, the Gauss-Seidel method is known to converge faster than the Jacobi method and, thus, it is preferable to use it. ⋆

This manuscript was originally written in July 2013, and the names and affiliation of the authors are reflecting the state at that time. Recently, the original manuscript went through a few rounds of AI-based reviews, which identified several typos and bugs, which were corrected, but which didn’t change the original ideas. AI was only used for identifying issues, not for correcting them.

However, when the size of the system of equations is large enough even sequential implementations of these methods are too slow and parallel versions are required. The Jacobi method can be easily parallelized, because the value of each unknown at the current iteration depends only on the values of the unknowns at the previous iteration. The Gauss-Seidel method, on the other hand, seems difficult to parallelize because the value of an unknown i at iteration k, x(k, i), depends on all the values x(k, j < i). This property is what helps the Gauss-Seidel method converge to a good solution with a smaller number of iterations than Jacobi, but it is also what prevents it from being easily parallelized. In this paper we propose a new parallel GaussSeidel algorithm for dense systems of linear equations which can overcome the limitations of the standard Gauss-Seidel method at the expense of more preprocessing time. In fact, with enough preprocessing time available and enough available processors, we can reach a parallel running time of O(log(n)) per iteration (n is the size of the system of linear equations), which is also the best result possible for the parallel version of the Jacobi method. However, this is not the main contribution of this paper. The most important result is that our algorithm can make use of a parameter m which can tune the amount of preprocessing time and the amount of parallelism available at each iteration of the Gauss-Seidel method. The algorithm will provide a trade-off between the preprocessing time and the best achievable parallel runtime: the more preprocessing time we can spend, the more parallelism per iteration we can achieve (and, thus, the better the parallel runtime can be). The rest of this paper is structured as follows. In Section 2 we describe the Gauss-Seidel method in detail and we discuss some parallelization opportunities of the standard formulation of this method. In Section 3 we present our new parallel Gauss-Seidel algorithm together with a detailed complexity analysis. The algorithm is focused on achieving the best parallel running time possible for a given number P of available processors. In Section 4 we conclude. 2. The Gauss-Seidel Method Gauss-Seidel is an iterative method for solving systems of linear equations A · x = b.

2

(1)

A is an n · n matrix, x and b are vectors with n elements, and x is the vector of unknowns. The method starts with an initial solution vector and after each iteration computes a new solution vector. The algorithm ends when the difference between the solution vector at the current iteration and the one from the previous iteration is below a predefined threshold. We define by x(0) the initial solution vector and by x(k) the solution vector obtained after k iterations. If we write A as A = D − L − U , where D is the diagonal part of A, −L is the strict lower part of A, and −U is the strict upper part of A, then we can write x(k + 1) as x(k + 1) = (D − L)−1 · U · x(k) + (D − L)−1 · b

(2)

If we denote by A(i, j) the element at row i and column j of A, by x(k, i) the element i of x(k) and by b(i) the element i of b (0 ≤ i, j ≤ n − 1), then we can compute x(k + 1) as: x(k + 1, i) = (b(i) −

X

A(i, j) · x(k + 1, j) −

j<i

X

A(i, j) · x(k, j))/A(i, i) (3)

j>i

The Gauss-Seidel method seems an inherently sequential method. x(k, i) can be computed only after all the values x(k, j < i) were already computed. The only obvious part which can be parallelized is the computation of each x(k, i) value (which consists of O(n) arithmetic operations). This computation (a reduction operation applied to n values) can make use of P ≤ n processors and can be implemented with a time complexity of O(n/P + log(P )) [10]. Thus, a Gauss-Seidel iteration can be very simply parallelized to run in O(n2 /P + n · log(P )) time but can use only at most n processors. For P = n the parallel runtime is O(n · log(n)). Note that in this paper we will use the terms and notations from [9] regarding the analysis of parallel algorithms (in particular, the terms parallel runtime and speedup). 3. The Parallel Gauss-Seidel Algorithm In [2] the authors were able to devise other more intelligent parallelization schemes than the obvious parallel algorithm mentioned in the previous section. These algorithms are suitable in GPU settings (where a large number of processors is available). However, none of the parallel algorithms presented in [2] has a parallel runtime better than O(n) per iteration no matter how 3

many processors are available. In this section we present a parallel GaussSeidel algorithm which is capable of achieving a sub-linear time complexity per iteration. However, we are capable of obtaining this time complexity at the expense of more preprocessing time. Note that obtaining the maximum degree of parallelism is not conceptually difficult. All we need to do is to compute E = (D − L)−1 · U and F = (D − L)−1 · b and then use Eq. 2. In Eq. 2 each component of x(k + 1) can be computed independently based only on the components of x(k). Thus, for P ≤ n processors we can obtain a parallel runtime per iteration of O(n2 /P ). For n < P ≤ n2 the achievable parallel runtime per iteration is O(n2 /P + log(P/n)). In order to obtain the parallel runtimes mentioned we need to first compute the inverse of the lower triangular matrix (D − L) and then multiply it by the upper triangular matrix U . Algorithms with O(n3 ) serial execution times exist for both operations. Parallel versions of these algorithms also exist. In [3] an algorithm with a parallel runtime of O(n3 /P ) for the inversion of triangular matrices has been proposed. However, the algorithm is only capable of using at most n processors for inverting an n · n matrix. An O(n) parallel runtime algorithm with O(n2 ) processors for inverting a triangular n · n matrix was obtained in [6]. The best parallel runtime that can be achieved by a known algorithm for inverting a triangular n · n matrix is O(log2 (n)) using O(n3 ) processors [4], [5]. The situation is better for parallel matrix multiplication algorithms [7]. There are several algorithms which can achieve a parallel runtime of O(n3 /P ) when the number of available processors P is at most n2 . An algorithm which can use n3 processors and achieve an O(log(n)) parallel running time was also analyzed in [7] and improved in order to achieve this running time with only n3 / log(n) processors. Actually, the mentioned algorithm has a parallel runtime of O(n3 /P + log(P/(n2 ))) when n2 ≤ P ≤ n3 processors are used. In general these computations may be prohibitive in terms of running time (in order to achieve very good parallel runtimes a large number of processors is required, which may not be available). The main contribution of this paper is to present a parallel Gauss-Seidel algorithm which can be tuned through a parameter m (1 ≤ m ≤ n). m controls the amount of parallelism available per iteration and the amount of preprocessing time - larger values of m lead to a larger degree of parallelism at the expense of more preprocessing time. In fact, m = 1 will denote a rather efficient Gauss-Seidel method with O(n) parallel runtime per iteration at best and m = n will represent the algorithm described earlier with the maximum degree of parallelism (O(log(n)) parallel 4

runtime per iteration at best). We will split the n unknowns into groups of size m. We will assume that m is a divisor of n. If not then we can increase n to the nearest multiple of m. We will also add extra rows and columns to the matrix A (all the elements on these extra rows and columns will be 0, except on the main diagonal where we will set the elements equal to 1) and extra elements to the vectors b and x. The original unknowns will not depend on the extra unknowns (whose values we are not interested in). We will denote by xg(k, i) the unknowns of group i after k iterations (xg(k, i) consists of the values x(k, m · i), . . . , x(k, m · (i + 1) − 1)). We define similarly bg(i) as the ith group of m consecutive values from the vector b. In the preprocessing stage we will compute an n/m by n/m matrix T , where each element is a block of size m·m. We will denote by T (i, j) the block at row i, column j. We will also introduce the vectors bgt(i), corresponding to the bg(i) vectors. With these notations we can rewrite Eq. 3 at block-level:

xg(k + 1, i) = bgt(i) +

X

T (i, j) · xg(k + 1, j) +

j<i

X

T (i, j) · xg(k, j)

(4)

j≥i

Basically, T (i, j) encodes the dependencies of xg(k + 1, i) on: (a) xg(k + 1, j) for j < i, or (b) xg(k, j) for j ≥ i. We will discuss the preprocessing stage and the implementation of the parallel Gauss-Seidel algorithm in two separate subsections. 3.1. The Preprocessing Stage An O(n2 · m) algorithm for computing all the entries in the matrix T can be designed as follows. For each group i(0 ≤ i ≤ n/m − 1), we will compute all the entries T (i, j)(0 ≤ j ≤ n/m − 1). First, we have T (i, i)(0, 0) = 0 and T (i, i)(0, 1 ≤ v ≤ m − 1) = −A(i · m, i · m + v)/A(i · m, i · m). Similarly, for j ̸= i, we have T (i, j)(0, 0 ≤ v ≤ m − 1) = −A(i · m, j · m + v)/A(i · m, i · m). Similarly, bgt(i)(0) = b(i · m)/A(i · m, i · m). For 1 ≤ u ≤ m − 1 we can compute T (i, j)(u, v) and bgt(i)(u) as follows: T (i, j)(j ̸= i, u ≥ 1, 0 ≤ v ≤ m − 1) = −A(i · m + u, j · m + v)/A(i · m + u, i · m + u)+ (5) Pu−1 w=0 (−A(i · m + u, i · m + w)/A(i · m + u, i · m + u)) · T (i, j)(w, v) 5

T (i, i)(u ≥ 1, 0 ≤ v ≤ u) = Pu−1 w=0 (−A(i · m + u, i · m + w)/A(i · m + u, i · m + u)) · T (i, i)(w, v)

(6)

T (i, i)(u ≥ 1, u < v ≤ m − 1) = −A(i · m + u, i · m + v)/A(i · m + u, i · m + u)+ Pu−1 w=0 (−A(i · m + u, i · m + w)/A(i · m + u, i · m + u)) · T (i, i)(w, v) (7) bgt(i)(u ≥ 1) = b(i · m + u)/A(i · m + u, i · m + u)+ Pu−1 w=0 (−A(i · m + u, i · m + w)/A(i · m + u, i · m + u)) · bgt(i)(w)

(8)

The computation of each value T (i, j)(u, v) and bgt(i)(u) can be easily parallelized to obtain a parallel runtime of O(m/P + log(P )) for P ≤ m processors. Since all the values of a single row can be computed independently, each row u of all the T (i, j) matrices (for fixed i and variable j) can be computed in O(m · n/P ) time if P ≤ n, or O(m · n/P + log(P/n)) time if n < P ≤ m·n. Since all the groups i(0 ≤ i ≤ n/m−1) are independent (there is no dependency between T (i, j) and T (i′ , j ′ ) with i ̸= i′ ), we can compute one row from each block in parallel (n/m rows in parallel), thus obtaining a parallel runtime of O(n2 /P ) if P ≤ n2 /m, or O(n2 /P + log(P · m/n2 )) if n2 /m < P ≤ n2 . In order to compute all the m rows for each group we need to multiply the time complexities mentioned earlier by a factor of m. A better approach when many processors P are available can be derived if we consider the (n/m) · (n/m) matrices A′ (i, j) of size m · m, where A′ (i, j)(u, v) = A(i · m + u, j · m + v) (0 ≤ u, v ≤ m − 1). If we let D′ (i), −L′ (i) and −U ′ (i) be the diagonal part, strict lower part and strict upper part, respectively, of A′ (i, i), then T (i, i) = (D′ (i) − L′ (i))−1 · U ′ (i) (this is similar to Eq. 2). Similarly, T (i, j) = (D′ (i)−L′ (i))−1 ·(−A′ (i, j)) (for j ̸= i), and bgt(i) = (D′ (i) − L′ (i))−1 · bg(i). In order to compute all the matrices T (i, j) (fixed i, 0 ≤ j ≤ n/m − 1) we need to compute the inverse of a lower triangular matrix and then multiply the result (also a lower triangular matrix) by an upper triangular matrix (U ′ (i)) or a square matrix (−A′ (i, j)). As mentioned earlier, there are sequential algorithms with O(m3 ) time complexities for both operations, but parallel algorithms also exist [3], [4], [7]. In the best case we can obtain an O(log2 (m)) parallel runtime with O(n2 · m) processors available for computing all the blocks of T . 6

3.2. The Implementation Once the matrix T is computed, we can use the following algorithm for the Gauss-Seidel method. We will compute the values of the unknowns at the current iteration in increasing order of their groups. For each group we can achieve the maximum degree of parallelism possible, but the computation of the values for a given group may only be finalized after all the values from the previous groups were fully computed. At each step of our algorithm we will need to perform O(n/m) independent matrix-vector multiplications (where each matrix has m rows and m columns and each vector has m elements). Parallel algorithms for the matrix-vector multiplication problem were considered in [8]. With P ≤ n processors available we can obtain a parallel runtime of O(n · m/P ) (for all the matrix-vector multiplications). With n < P ≤ n · m processors available we can obtain a parallel runtime of O(n·m/P +log(P/n)). Each iteration of our Gauss-Seidel algorithm consists of n/m steps. Thus, the parallel runtime per iteration is O(n2 /P ) for P ≤ n, or O(n2 /P + (n/m) · log(P/n)) for n < P ≤ n · m. Besides the parallel matrix-vector multiplications there are several other simpler operations which need to be performed by the algorithm (e.g. initializing a vector in parallel, copying a vector in parallel, computing the difference of two vectors, computing the sum of the elements of a vector, etc.), which have parallel runtimes which are better than the parallel matrixvector multiplication for the same number of available processors. Thus, the time required to perform matrix-vector multiplications dominates the running time of our algorithm. An iteration of our parallel Gauss-Seidel method is presented in Algorithm 1 while the full algorithm is presented in Algorithm 2. The algorithm maintains for each group i two vectors of m elements: xgcurr(i) and xgold(i). These vectors obey the following invariants. At the beginning of the step i of the current iteration (k + 1) we have: xgcurr(i) already contains the values corresponding to xg(k + 1, i) from Eq. 4, and xgold(i) contains the values of the unknowns from the previous iteration (xg(k, i) from Eq. 4). At the end of the step i, xgold(i) ends up containing the values of the unknowns at the end of the iteration (xgcurr(i) is copied over into xgold(i)). Note that Algorithm 1 follows a slightly different formulation than the one in Eq. 4, because when step i starts, the values xg(k + 1, i) are already computed, while Eq. 4 naturally implies that these values would need to be computed during the execution of the step.

7

Algorithm 1 Parallel Gauss-Seidel Iteration ε=0 for i = 0 to n/m − 1 do Do in parallel delta(i) P= xgcurr(i) − xgold(i) Do in parallel ε += j=0..m−1 |delta(i)(j)| Do in parallel xgold(i) = xgcurr(i) for j = 0 to n/m−1 do {Each iteration j can be executed independently, in parallel.} Compute in parallel tmp(j, i) = T (j, i) · delta(i) Do in parallel xgcurr(j) += tmp(j, i) end for end for

Algorithm 2 Parallel Gauss-Seidel Algorithm Compute in parallel the matrix T and the vectors bgt(i) (0 ≤ i ≤ n/m − 1) Do in parallel xgold(i) = 0 (0 ≤ i ≤ n/m − 1) Do in parallel xgcurr(i) = bgt(i) (0 ≤ i ≤ n/m − 1) ε=∞ while ε > εconvergence do Parallel Gauss-Seidel Iteration end while

8

4. Conclusions In this paper we presented a novel parallel Gauss-Seidel algorithm for solving dense systems of linear equations. The algorithm can obtain better parallel execution times than previous algorithms at the expense of more preprocessing time. The algorithm’s preprocessing time can be tuned, so that the more preprocessing time we can accept, the more parallelism per iteration we can obtain. Our algorithm may also be of particular interest when we need to solve multiple systems of linear equations which have the same matrix A but different vectors b. Computing the inverses (D′ (i) − L′ (i))−1 for each block A′ (i, i) and the matrix T is the most time-consuming part of the preprocessing stage and only depends on the matrix A, so it only needs to be performed once. For each new vector b, we only need to compute the vectors bgt(i) = (D′ (i) − L′ (i))−1 · bg(i) via matrix-vector multiplications (or via Eq. 8). References [1] J. D. Faires, R. L. Burden, “Numerical Methods”, 4th edition, Brooks/Cole Cengage Learning, 2013. [2] H. Courtecuisse, J. Allard, “Parallel Dense Gauss-Seidel Algorithm on Many-Core Processors”, Proceedings of the High Performance Computation Conference (HPCC), pp. 139-147, 2009. [3] W. Nasri, Z. Mahjoub, “Design and Implementation of a General Parallel Divide and Conquer Algorithm for Triangular Matrix Inversion”, International Journal of Parallel and Distributed Systems and Networks, Vol. 5(1), 2002. [4] A. H. Sameh, R. P. Brent, “Solving Triangular Systems on a Parallel Computer”, SIAM J. Numer. Anal., Vol. 14 (6), pp. 1101-1113, 1977. [5] A. Schikarski, D. Wagner, “Efficient Parallel Matrix Inversion on Interconnection Networks”, Journal of Parallel and Distributed Computing, Vol. 34 (2), pp. 196-201, 1996. [6] A. Kiper, “An Efficient Parallel Triangular Inversion by Gauss Elimination with Sweeping”, Proceedings of the 4th International Euro-Par Conference on Parallel Processing, pp. 793-797, 1998. 9

[7] A. Gupta, V. Kumar, “Scalability of Parallel Algorithms for Matrix Multiplication”, Proceedings of the International Conference on Parallel Processing (ICPP), pp. 115-123, 1993. [8] A. Grama, A. Gupta, G. Karypis, V. Kumar, “Introduction to Parallel Computing”, 2nd edition, Pearson Education, 2003. [9] T. Rauber, G. Rünger, “Parallel Programming for Multicore and Cluster Systems”, Springer, 2010. [10] D. Skillicorn, “Foundations of Parallel Programming”, Cambridge University Press, 1994.

10

Related documents

• An Exact Parallel Dense Gauss-Seidel Algorithm with Trade-Offs between Parallel Execution Time and Preprocessing Time#1491598
• Exploiting Task-Based Parallelism for the Red-Black Gauss-Seidel Method on 2D Grids#332495
• How the Dirty Set Moves: Characterizing Schedule Reuse in Incremental Sparse Matrix Multiplication#1345312
• How the Dirty Set Moves: Characterizing Schedule Reuse in Incremental Sparse Matrix Multiplication#1335967
• Asynchronous Jacobi and randomized Gauss--Seidel methods in shared and distributed memory: A unified convergence-rate analysis#1109591