Thomas algorithm, pointwise version

Thomas algorithm, pointwise version
Sequential algorithm
Serial complexity	[math]8n-7[/math]
Input data	[math]4n-2[/math]
Output data	[math]n[/math]
Parallel algorithm
Parallel form height	[math]5n-4[/math]
Parallel form width	[math]2[/math]

Primary authors of this description: A.V.Frolov, Vad.V.Voevodin (Section 2.2), Alexey Teplov (Sections 2.4, 2.5)

1 Properties and structure of the algorithm

1.1 General description of the algorithm

Thomas algorithm is a variant of Gaussian elimination used for solving a system of linear algebraic equations (SLAE) [math]Ax = b[/math], where

[math] A = \begin{bmatrix} a_{11} & a_{12} & 0 & \cdots & \cdots & 0 \\ a_{21} & a_{22} & a_{23}& \cdots & \cdots & 0 \\ 0 & a_{32} & a_{33} & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & a_{n-1 n-2} & a_{n-1 n-1} & a_{n-1 n} \\ 0 & \cdots & \cdots & 0 & a_{n n-1} & a_{n n} \\ \end{bmatrix}, x = \begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \\ \end{bmatrix}, b = \begin{bmatrix} b_{1} \\ b_{2} \\ \vdots \\ b_{n} \\ \end{bmatrix} [/math]

However, presentations of the elimination method often use a different notation and a numbering for the right-hand side and matrix of the system. For instance, the above SLAE can be written as

[math] A = \begin{bmatrix} c_{0} & -b_{0} & 0 & \cdots & \cdots & 0 \\ -a_{1} & c_{1} & -b_{1} & \cdots & \cdots & 0 \\ 0 & -a_{2} & c_{2} & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & -a_{N-1} & c_{N-1} & -b_{N-1} \\ 0 & \cdots & \cdots & 0 & -a_{N} & c_{N} \\ \end{bmatrix}\begin{bmatrix} y_{0} \\ y_{1} \\ \vdots \\ y_{N} \\ \end{bmatrix} = \begin{bmatrix} f_{0} \\ f_{1} \\ \vdots \\ f_{N} \\ \end{bmatrix} [/math]

(here, N=n-1). If each equation is written separately, then we have

[math]c_{0} y_{0} - b_{0} y_{1} = f_{0}[/math],

[math]-a_{i} y_{i-1} + c_{i} y_{i} - b_{i} y_{i+1} = f_{i}, 1 \le i \le N-1[/math],

[math]-a_{N} y_{N-1} + c_{N} y_{N} = f_{N}[/math].

Here, we examine the so-called right-elimination, which is the variant of elimination method in which a SLAE is processed top-down and then in the reverse direction. In this variant, one first goes top down eliminating the subdiagonal unknowns and then bottom up eliminating the superdiagonal ones. There is also a variant, called the left-elimination, in which a SLAE is first processed bottom-up and then top-down. There is basically no difference between both variants; consequently, we do not include a separate description of the left-elimination.

Figure 1. Graph of the elimination method for n=8. The input and output data are not shown. The symbols / and f denote division and operation a+bc or a-bc, respectively

1.2 Mathematical description of the algorithm

In the notation introduced above, the forward elimination path consists in calculating the elimination coefficients

[math]\alpha_{1} = b_{0}/c_{0}[/math],

[math]\beta_{1} = f_{0}/c_{0}[/math],

[math]\alpha_{i+1} = b_{i}/(c_{i}-a_{i}\alpha_{i})[/math], [math]\quad i = 1, 2, \cdots , N-1[/math],

[math]\beta_{i+1} = (f_{i}+a_{i}\beta_{i})/(c_{i}-a_{i}\alpha_{i})[/math], [math]\quad i = 1, 2, \cdots , N[/math].

Then, in the backward elimination path, one calculates the solution

[math]y_{N} = \beta_{N+1}[/math],

[math]y_{i} = \alpha_{i+1} y_{i+1} + \beta_{i+1}[/math], [math]\quad i = N-1, N-2, \cdots , 1, 0[/math].

According to the literature (see ), these formulas are equivalent to calculating a variant of the [math]LU[/math] decomposition for the coefficient matrix accompanied by solving bidiagonal systems via forward and back substitutions.

1.3 Computational kernel of the algorithm

The computational kernel of this algorithm can be thought of as compiled of two parts, namely, the forward and backward elimination paths. The computational kernel of the forward elimination path consists of sequences of divisions, multiplications, and additions/subtractions. The computational kernel of the backward elimination path contains only multiplications and additions sequences.

Figure 2. Detailed graph of the elimination method for n=8 (each reciprocal number is only once calculated). The input and output data are not shown. Each calculation of a reciprocal number is marked by inv, and each multiplication is marked by mult. The operations repeated after the replacement of the right-hand side of SLAE are set out in grey

1.4 Macro structure of the algorithm

The macro structure of this algorithm can be represented as the combination of the forward and backward elimination paths. In addition, the forward path can also be split into two macro units, namely, the triangular decomposition of the coefficient matrix and the forward substitution for a bidiagonal SLAE. These units are executed simultaneously, that is, in parallel. The process of solving the bidiagonal SLAE uses the results produced by the triangular decomposition.

1.5 Implementation scheme of the serial algorithm

The method is executed as the following sequence of steps:

1. Initialize the forward elimination path:

[math]\alpha_{1} = b_{0}/c_{0}[/math],

[math]\beta_{1} = f_{0}/c_{0}[/math]

2. For [math]i[/math] increasing from [math]1[/math] to [math]N-1[/math], execute the forward path formulas:

[math]\alpha_{i+1} = b_{i}/(c_{i}-a_{i}\alpha_{i})[/math],

[math]\beta_{i+1} = (f_{i}+a_{i}\beta_{i})/(c_{i}-a_{i}\alpha_{i})[/math].

3. Initialize the backward elimination path:

[math]y_{N} = (f_{N}+a_{N}\beta_{N})/(c_{N}-a_{N}\alpha_{N})[/math]

4. For [math]i[/math] decreasing from [math]N-1[/math] to [math]0[/math], execute the backward path formulas: [math]y_{i} = \alpha_{i+1} y_{i+1} + \beta_{i+1}[/math].

The formulas of the forward path contain double divisions by the same expressions. These divisions can be replaced by the calculation of reciprocal numbers succeeded by the multiplication by these numbers.

1.6 Serial complexity of the algorithm

Consider a tridiagonal SLAE consisting of [math]n[/math] equations with [math]n[/math] unknowns. The elimination method as applied to solving such a SLAE in its (fastest) serial form requires

• [math]2n-1[/math] divisions,
• [math]3n-3[/math] additions/subtractions, and
• [math]3n-3[/math] multiplications.

Thus, in terms of serial complexity, the elimination method is qualified as a linear complexity algorithm.

1.7 Information graph

The information graph of the elimination method is shown in fig. 1. An analysis indicates that this graph is virtually serial. During the forward path, two branches (namely, the left branch, corresponding to the matrix decomposition, and the central branch, responsible for the solution of the first bidiagonal system) can be executed in parallel. The right branch corresponds to the backward path. It is evident from the figure that not only the mathematical content of this process but also the structure of the algorithm graph and the direction of the corresponding data flows are in complete agreement with the appellation "backward path." The variant with divisions replaced by reciprocal number calculations is illustrated by the graph in fig. 2.

1.8 Parallelization resource of the algorithm

The parallel version of the elimination method as applied to solving a tridiagonal SLAE consisting of [math]n[/math] equations with [math]n[/math] unknowns requires that the following layers be executed:

• [math]n[/math] division layers (all the layers except for one contain [math]2[/math] divisions),
• [math]2n - 2[/math] multiplication layers and the same number of addition/subtraction layers ([math]n-1[/math] layers contain two operations of both types, while another [math]n-1[/math] layers contain a single operation of each type).

Thus, in terms of the parallel form height, the elimination method is qualified as an algorithm of complexity [math]O(n)[/math]. In terms of the parallel form width, its complexity is [math]2[/math].

1.9 Input and output data of the algorithm

Input data: tridiagonal matrix [math]A[/math] (with entries [math]a_{ij}[/math]), vector [math]b[/math] (with components [math]b_{i}[/math]).

Output data: vector [math]x[/math] (with components [math]x_{i}[/math]).

Size of the output data: [math]n[/math].

1.10 Properties of the algorithm

It is clearly seen that the ratio of the serial to parallel complexity is a constant (which is less than [math]2[/math]).

The computational power of the algorithm, understood as the ratio of the number of operations to the total size of the input and output data, is also a constant.

Within the framework of the chosen version, the algorithm is completely determined.

Routinely, the elimination method is used for solving SLAEs with diagonally dominant coefficient matrices. For such systems, the algorithm is guaranteed to be stable.

Suppose that several systems with the same coefficient matrix have to be solved. Then the left branch of calculations (see the figures with the algorithm graph) may not to be repeated. The reason is that the [math]LU[/math] decomposition of the coefficient matrix needs not to be recalculated. In such cases, the variant with divisions replaced by reciprocal number calculations is preferable.

2 Software implementation of the algorithm

2.1 Implementation peculiarities of the serial algorithm

Different implementations of the algorithm are possible depending on how the coefficient matrix is stored (as a single array with three rows or three distinct arrays) and how the calculated coefficients are stored (in locations of the already used matrix entries or separately).

Let us give an example of a Fortran subroutine implementing the elimination method. Here, all the matrix entries are stored in a single array; moreover, the entries neighboring in a row are located closely to each other, and the calculated coefficients are stored in locations of the original matrix entries.

      subroutine progm (a,x,N)
      real a(3,0:N), x(0:N)
      a(2,0)=1./a(2,0)
      a(3,0)=-a(3,0)*a(2,0) ! alpha 1
      x(0)=x(0)*a(2,0) ! beta 1
      do 10 i=1,N-1
        a(2,i)=1./(a(2,i)+a(1,i)*a(2,i-1))
        a(3,i)=-a(3,i)*a(2,i) ! alpha i+1
        x(i)=(x(i)-a(1,i)*x(i-1))*a(2,i) ! beta i+1
  10  continue
      a(2,N)=1./(a(2,N)+a(1,N)*a(2,N-1))
      x(N)=(x(N)-a(1,N)*x(N-1))*a(2,N) ! y N
      do 20 i=N-1,0,-1
        x(i)=a(3,i)*x(i+1)+x(i) ! y i
  20  continue
      return
      end

There are many simple implementations of the method on the web; see, for instance, Wikibooks - Algorithm implementation - Linear Algebra - Tridiagonal matrix algorithm.

2.2 Locality of data and computations

One can see from the algorithm graph that the spatial locality of the data is fair. Indeed, all the arguments needed for an operation are calculated "nearby." However, the temporal locality is not as good. Suppose that the cache cannot hold all the problem data. Then repeated cache misses will occur while performing calculations in the upper left corner of a SLAE. Thus, one of the recommendations to application engineers using the elimination method is as follows: Arrange the calculations so that all the elimination paths are sufficiently short and their data can be placed in the cache.

The elimination method is an example of a serial algorithm with the locality of computations so high that it even becomes excessive . Computational cores of processors virtually cannot use their superscalar capabilities because the data needed by the basic operations of this method are produced by immediately preceding operations. This sharply deteriorates the performance of the elimination method even for modern single-processor and single-core systems.

2.2.1 Locality of implementation

2.2.1.1 Structure of memory access and a qualitative estimation of locality

Figure 3. Elimination method, pointwise version. The general memory access profile

Figure 3 presents the memory access profile for an implementation of the pointwise version of elimination method. This profile is formed of accesses to 5 arrays storing three diagonals of the coefficient matrix, the vector of right-hand sides, and the resulting solution. Judging from the general profile, one can say that the subroutine is composed of two stages. At the first stage, it performs the serial access to entries of four arrays. The second stage is also the serial access performed in the reverse order.

Figure 4. An isolated fragment of the general profile (set out in green in fig. 3)

However, a closer look shows that the structure of the profile is somewhat more complicated (see fig. 4). At each step, some of the arrays are accessed several times to very close locations, which improves both spatial and temporal locality. Thus, the general profile has a high spatial locality (because serial accesses of the entries dominate here) and a medium temporal locality (because certain entries are immediately used again).

2.2.1.2 Quantitative estimation of locality

The basic fragment of the implementation used for obtaining quantitative estimates is given here (Kernel function). The start-up parameters are described here.

The first estimate is produced on the basis of daps, which estimates the number of memory accesses (reading and writing) per second. This characteristic, used for the interaction with memory, is an analog of the flops estimate. It largely estimates the performance of this interaction rather than locality. Nevertheless, this characteristic (in particular, its comparison with the next characteristic cvg) is a good source of information.

Figure 5 presents the values of daps for implementations of popular algorithms. They are arranged in increasing order (in general, the larger daps, the higher efficiency). The daps value of the subroutine under discussion is fairly unexpected. According to the above qualitative analysis, the locality of memory access is sufficiently high for this subroutine. However, daps estimates the performance of memory interaction to be worse than for the RandomAccess test! Such a distinction between locality and performance occurs not very often.

A detailed inspection of the original code reveals two reasons of this phenomenon. First of all, the iterations of the main loop are memory interdependent; that is, each iteration uses the data calculated by the preceding iteration. In this case, such an interdependence seems to significantly deteriorate the performance of memory interaction; on the other hand, it is rather good from the viewpoint of locality.

The second discovery is that, unexpectedly, one of divisions considerably slows down the execution of the subroutine. After this division was replaced by multiplication, the execution time reduced by a factor of 2.5. The assembler codes of two versions looked identical (except for the obvious replacement of a single division by multiplication). These features of the subroutine are not reflected in its locality properties, which resulted in such a substantial gap between the locality estimate and the performance of memory interaction.

Figure 5. Comparison of the values of daps

The second characteristic, called cvg, is intended for obtaining a locality estimate that would be more machine independent. It determines how often a program needs to fetch data to the cache memory. Accordingly, smaller values of cvg correspond to less often fetch operations and better locality.

Figure 6 presents the values of cvg for the same set of implementations. They are arranged in decreasing order (in general, the smaller cvg, the higher locality). Here, all our previous conclusions are confirmed by cvg: Locality turns out to be fairly good, which was already noted when discussing the qualitative locality estimate.

Figure 6. Comparison of the values of cvg

2.3 Possible methods and considerations for parallel implementation of the algorithm

As seen from the analysis of the algorithm graph, it virtually cannot be parallelized (if not substantially modified). Consequently, there are two ways of adapting eliminations for parallel computers. If a problem uses eliminations, then the first way is to split it so as to obtain fairly many of them and assign a separate processor (separate core) to each elimination. The second way is to replace the elimination method by its parallel variants (the cyclic reduction method, serial-parallel variants, etc.).

2.4 Scalability of the algorithm and its implementations

There is no sense in discussing the scalability of the elimination method because it is completely non-parallel. However, we note that, for its parallel variants, the scalability should be analyzed in comparison with a single processor implementation of the classical method rather than with single processor implementations of these parallel variants.

An analysis of scalability was performed for the algorithm implemented in accordance with the methodology proposed in http://algowiki-project.org/ru/Scalability_methodology. The analysis was based on the experiments conducted on the Lomonosov supercomputer of the Moscow University Supercomputing Center].

The variable start-up parameters of this implementation and their ranges are as follows:

• the number of processors is 1;
• the domain size varies in the range [10240 : 1024000] with the step 10240.

The experiments resulted in the following range for the efficiency of this implementation:

• the minimum efficiency is 0.019%;
• the maximum efficiency is 0.0255%.

The following two figures show how the performance and efficiency of this implementation depend on the variable start-up parameters.

Figure 7. Implementation of the algorithm. Performance as a function of the vector size

Figure 8. Implementation of the algorithm. Efficiency as a function of the vector size

The low efficiency is likely to be related to the excessive locality discussed in the section on the locality of data and calculations.

The implementation discussed in this section is presented here.

2.5 Dynamic characteristics and efficiency of the algorithm implementation

In view of the serial nature of the algorithm and its excessive locality, the analysis of its dynamic characteristics seems to be of little value; therefore, it was not conducted.

2.6 Conclusions for different classes of computer architecture

The elimination method is intended for the classical von Neumann architecture. It is not even suitable for the effective loading of single core systems supporting superscalarity. In order to parallelize the process of solving a SLAE with a tridiagonal coefficient matrix, one should use some parallel substitute of this method, for instance, the most popular cyclic reduction or the new serial-parallel version of the elimination method. The latter yields to the former in terms of the critical path of the graph, but, instead, its graph has a more regular structure.

2.7 Existing implementations of the algorithm

The elimination method is present in all the standard program packages for linear algebra problems. However, it is so simple that most of the users (researchers and application engineers) prefer to write their own program segments implementing this method. It should be noted, though, that the elimination method in the standard packages is ordinarily not implemented in its pure form but rather as the pair "factorization into the product of two bidiagonal matrices" and "solution of a SLAE with the already factorized tridiagonal coefficient matrix." This is, for instance, the case with the SCALAPACK package and its predecessors.