A Novel Parallel Algorithm Based on the Gram-Schmidt Method for Tridiagonal Linear Systems of Equations

Abstract

This paper introduces a new parallel algorithm based on the Gram-Schmidt orthogonalization method. This parallel algorithm can find almost exact solutions of tridiagonal linear systems of equations in an efficient way. The system of equations is partitioned proportional to number of processors, and each partition is solved by a processor with a minimum request from the other partitions' data. The considerable reduction in data communication between processors causes interesting speedup. The relationships between partitions approximately disappear if some columns are switched. Hence, the speed of computation increases, and the computational cost decreases. Consequently, obtained results show that the suggested algorithm is considerably scalable. In addition, this method of partitioning can significantly decrease the computational cost on a single processor and make it possible to solve greater systems of equations. To evaluate the performance of the parallel algorithm, speedup and efficiency are presented. The results reveal that the proposed algorithm is practical and efficient.

1. Introduction

Linear systems of equations occur frequently in many fields of engineering and mathematics. Different approaches have been used to solve the linear systems of equations [1], for example, 𝑄𝑅 factorization. There are different 𝑄𝑅 factorization methods, such as Givens transformations, Householder transformations, and Gram-Schmidt orthogonalisation.

A considerable amount of study has been done on parallel 𝑄𝑅 factorization algorithms [2, 3]. In addition, some numerical libraries and softwares have supplied the 𝑄𝑅 factorization, such as ScaLAPACK [4], and Parallel Linear Algebra Software for Multicore Architectures (PLASMA) [5]. PLASMA Tile 𝑄𝑅 factorization is provided for multicore system architectures using the block algorithms [6].

Tridiagonal system is a well known form that arises in areas such as linear algebra, partial differential equations, orthogonal polynomials, and signal processing. Problem of solving tridiagonal systems on parallel computers has been studied widely. It seems that first parallel tridiagonal system solver, referred to as cyclic reduction, was introduced by Hockney in 1965 [7]. Throughout more than forty years since then, different parallel algorithms have been proposed and developed. The recursive doubling algorithm was presented by Stone for fine grained parallel systems in 1973 [8]. Wang suggested a partitioning algorithm for coarse grained systems in 1981 [9]. The classical method is the inverse of the matrix containing unknown elements; recently, an explicit formula was given by Kilic [10].

In this paper, we consider a novel parallel algorithm based on the modified Gram-Schmidt orthogonalization method [11]. This method is a kind of 𝑄𝑅 factorization, in which firstly orthogonal and upper-triangular matrices are computed and then the solution of the system of equations is found using these matrices. Amodio et al. proposed a parallel factorization for tridiagonal matrices [12, 13]. On the other hand, parallel implementation of modified Gram-Schmidt algorithm has been studied in different literatures [14–17]. Although the proposed method in this paper is currently developed and verified for tridiagonal systems, the authors are working on extension of the method to other systems.

The scalability of parallel algorithms is a vital problem, so different techniques have been developed to overcome it such as processors load balancing [18] and reducing the number of synchronizations [19]. The obtained results show that the suggested parallel algorithm is practically scalable.

Therefore, the objective of this paper is to propose an efficient parallel algorithm. In Section 2, the 𝑄𝑅 factorization is introduced as a base of this algorithm. In Section 3, the procedure of the suggested parallel algorithm is completely described. A numerical example is given in Section 4. In Section 5, accuracy, speedup, and efficiency of the algorithm are evaluated. Finally, conclusions are given in Section 6.

2. 𝑄𝑅 Factorization

A 𝑄𝑅 factorization of a nonsingular matrix 𝐴∈ℝ_𝑛×𝑛is a decomposition of the matrix into a unitary orthogonal matrix 𝑄∈ℝ_𝑛×𝑛 (i.e.,𝑄^𝑇𝑄=𝐼) and an upper-triangular matrix 𝑅∈ℝ^𝑛×𝑛as𝐴=𝑄𝑅, (2.1) where the decomposition is unique. One of the most well-known 𝑄𝑅 methods is the Classical Gram-Schmidt (CGS) process. To overcome numerical instability, the modified Gram-Schmidt algorithm (MGS) is preferred after a simple change in the classical algorithm [20]. Each column of orthogonal matrix can be normalized as follows:𝑞_𝑘=𝑢^(𝑘)𝑢^(𝑘), (2.2) where 𝑞^(𝑘) is a column of orthogonal matrix and 1≤𝑘≤𝑛. Column 𝑢^(𝑘) is determined from previously computed orthogonal columns 𝑞^(𝑖), where 1≤𝑖≤𝑘−1, as𝑢^(𝑘)=𝑎^(𝑘)−^𝑘−1_𝑖=1𝑎^(𝑘),𝑞^(𝑖)𝑞^(𝑖), (2.3) where (𝑎^(𝑘),𝑞^(𝑖))𝑞^(𝑖) denotes the projection operator of the column 𝑎^(𝑘) onto the column 𝑞^(𝑖) and 𝑎^(𝑘),𝑞^(𝑖)=𝑎^(𝑘),𝑞^(𝑖)𝑞^(𝑖),𝑞^(𝑖), (2.4) where, ⟨𝑎^(𝑘),𝑞^(𝑖)⟩ presents the inner product of the columns𝑎^(𝑘) and 𝑞^(𝑖). Loss of orthogonality depends on matrix 𝐴, and the amount of deviation from orthogonality is minimized when 𝐴 is a well-conditioned matrix [21].

Additionally, an upper-triangular matrix is constructed column by column as𝑟_𝑖𝑘=𝑢^(𝑘),𝑞^(𝑖), (2.5) where 𝑟_𝑖𝑘 is a typical element of upper-triangular matrix and 1≤𝑖≤𝑛.

Furthermore, the 𝑄𝑅 factorization can be used as a part of solution of a linear system of equations. To obtain the solution of a linear system𝐴𝑋=𝐵,[7]𝑋∈ℝ^𝑛,𝐵∈ℝ^𝑛, (2.6) we consider the 𝑄𝑅 factorization 𝑄𝑅𝑋=𝐵. (2.7) Then we multiply the transpose of 𝑄 by both sides of (2.7), where 𝑄^𝑇𝑄=𝐼 and 𝑄^𝑇𝐵=𝐵^′. Finally, a linear system with an upper-triangular coefficient matrix is achieved as𝑅^′𝑋=𝐵^′, (2.8) where 𝑅^′=𝑅. The final step is computing the matrix 𝑋 by usual back substitution.

3. Parallel Algorithm

The concept of novel parallel algorithm based on the Gram-Schmidt method will be described here. The parallel implementation of the Gram-Schmidt method is not efficient compared with other methods, particularly when it is employed to solve a tridiagonal system of equations. However, some modifications allow appropriate conditions needed to parallelize the Gram-Schmidt method. In the Gram-Schmidt method, computation of each column of the orthogonal and upper-triangular matrices depends on the previously computed columns of the orthogonal matrix.

Most of the elements of a tridiagonal matrix are zero, except for those on three diagonals. Each term of projection operator in (2.3) consists of inner product of one column of this tridiagonal matrix with a computed column of orthogonal matrix. When most of the elements of one of these columns are zero, then a few elements of projection operator can be nonzero. As a result, only the last two terms of projection operator in the computation of each column of orthogonal matrix are nonzero, as 𝑟_𝑖𝑘≠0,𝑖=𝑘−2,𝑘−1. (3.1)

Figure 1 shows the terms of (2.3) for the first six columns of orthogonal matrix𝑄.

Figure 1 

The zero and nonzero terms in the equation of each orthogonal column, where 𝑞(𝑚) is the 𝑚th column of orthogonal matrix and 𝑎(𝑚) is 𝑚th column of matrix 𝐴. Also, the zero terms are crossed out.

As a result, computation of each column of the orthogonal and upper-triangular matrices depends only on the two previous adjacent columns of the orthogonal matrix𝑄. It means that, if the column (𝑘) of orthogonal matrix is computing, only two computed columns (𝑘−1) and (𝑘−2) of matrix 𝑄 and the column (𝑘) of matrix 𝐴 are requested, as𝑢^(𝑘)=𝑎^(𝑘)−𝑎^(𝑘),𝑞^(𝑘−2)𝑞^(𝑘−1)−𝑎^(𝑘),𝑞^(𝑘−1)𝑞^(𝑘−1). (3.2) This property is a base of our parallel algorithm. In this situation, if we omit two adjacent columns of matrix 𝐴, then all terms of projection operator for the next columns are zero. It means that the relation of columns of orthogonal matrix in (2.3) is disrupted. In order to use this feature, if two columns of matrix𝐴_𝑛×𝑛, for instance, 𝑎^(𝑘−2) and𝑎^(𝑘−1), are transferred to the end of matrix, the next corresponding columns in the orthogonal matrix are disconnected from the previous ones. It means that (𝑎^(𝑘),𝑞^(𝑘−2))=0 and(𝑎^(𝑘),𝑞^(𝑘−1))=0 and as a result,𝑢^(𝑘)=𝑎^(𝑘). (3.3) This transformation decreases the amount of computation significantly.

The parallel algorithm is achieved by partitioning the matrix 𝐴_𝑛×𝑛 into 𝑚−1 (𝑚 is the number of processors) blocks of fairly equal size, that is, each block contains approximately 𝑁_𝑒=𝑛/(𝑚−1) columns, and each processor is concerned with computing only its block. The computation can be done mostly independently by each processor, but at some points, communication is required. The first column of the matrix owned by processor 𝑒 is denoted as𝑗_𝑒, where 𝑗_𝑒<𝑗_𝑒+1 for𝑒=1,…,𝑚−1 and𝑁_𝑒=𝑗_𝑒+1−𝑗_𝑒.

To disconnect the relationship between adjacent blocks, two columns between them, that is, 𝑗_𝑒+1−1 and𝑗_𝑒+1 for𝑒=1,…,𝑚−2, are transferred to the end of matrix (3.4). These transferred columns make another partition, as shown in Figure 2. Therefore, the number of columns of the new block is𝐻=2×(𝑚−2): 𝑗_𝑚+(𝑒−1)⟵𝑗_𝑒+1−1,𝑗_𝑚+𝑒⟵𝑗_𝑒+1for𝑒=1,…,𝑚−2. (3.4) It should be noticed that the place of each element of matrix𝑋 corresponding to the transferred columns of matrix 𝐴 is also changed.

(a)

(b)

Figure 2 

The partitioning (a) and transformation (b) of matrix's columns. The amount of partition can be related to the number of processor. The last partition is made by transferred columns from between other blocks.

Each processor, except for the last one, can start computing some columns of the orthogonal and upper-triangular matrices that are related to the columns of its own block without receiving any data from other processors. This feature plays a significant role in improving the speedup of the algorithm. If processor 𝑃_𝑒 is responsible for a block consisting of columns from 𝑎^(𝑗𝑒) to𝑎^{(𝑗𝑒+1−1)}, it computes the columns from 𝑞^(𝑗𝑒) to 𝑞^{(𝑗𝑒+1−1)} of the orthogonal matrix with a minimum request for data of other partition, and the same holds for 𝑟^(𝑗𝑒) to 𝑟^{(𝑗𝑒+1−1)} of the upper-triangular matrix. All the blocks, except the last one, only need data of the last block, because the first two and last two columns of each block have some relation to some columns of the last block.

To achieve more speedup, another switching is used in each block, except for the last block. The last two columns in each block are transferred between the second and third columns of block. This switching causes the first four columns to be computed firstly and then be sent to the last processor before other columns of block are determined. As a result, all processors are in use simultaneously, and the efficiency of the parallel algorithm increases.

Lastly, all columns of the orthogonal and upper-triangular matrices are computed with the least network transaction. To show the final result, each processor sends its block to the last processor.

The next step is to find the solution of the linear system of equations. According to (2.7), the values of 𝑄^𝑇𝑄 and 𝑄^𝑇𝐵 should be determined. Because data of each block of matrix𝑄 and also matrix𝐵 exist in each processor, the computation of these matrices can be done without any request to the other processors' data.

Finally, the unknown matrix is obtained by a back substitution procedure. In this step, each processor just needs to communicate with the last processor. The back substitution is started from the last processor which consists of 𝐻 columns. The last 𝐻 elements of the unknown matrix are computed using just the last processor and are sent to other ones. Because the 𝐻 elements of matrix 𝑋 are computed, the last processor can compute independently, the last 𝐻 terms of each linear equation. As a result, instead of 𝐻 columns, just one column is needed to broadcast to the other processors. Then the other processors compute the related unknown elements with the received data. The complete flowchart of described parallel algorithm is shown in Figure 3.

Figure 3 

The flowchart of proposed parallel algorithm.

It may seem that the amount of computation in the last processor is more than the others. However, it should be mentioned that the number of columns in the last block is less than the number of columns in the others. If the number of unknowns increases, the difference between size of the last block and other blocks becomes greater.

4. Example

Let 𝐴 be a 15 × 15 nonsingular tridiagonal matrix. The results of three cases are presented here and compared with each other; we use the label case (a) for one processor, case (b) for two, and case (c) for three. In Figure 4, the transformations of two columns in case (b) (two processors) and four columns in case (c) (three processors) are shown.

(a)

(b)

(c)

Figure 4 

The matrix𝐴15×15; (a) one processor, (b) two processors, and (c) three processors.

The computed orthogonal matrices in the three cases are presented in Figure 5. The number of nonzero elements is 135 in case (a), 97 in case (b), and “93” in case (c). Hence, the computational cost decreases as the computational power (i.e., the number of processors) increases. This situation is an ideal condition for a parallel algorithm.

(a)

(b)

(c)

Figure 5 

The orthogonal matrices; (a) one processor, (b) two processors, and (c) three processors.

Finally, the upper-triangular matrices are shown in Figure 6. Notice that, in cases (b) and (c), the elements in the last row of each block are zero, except for elements that belong to the last processors, in addition to an element which is equal to 1. Therefore, the back substitution procedure in each block does not require the data of the other processors, only the last processor. This feature is also significant for solving a huge system of equations.

(a)

(b)

(c)

Figure 6 

The Upper-triangular matrices; (a) one processor, (b) two processors, and (c) three processors.

5. Results

The testbed specification is shown in Table 1. The results of using up to 16 processors are presented here. The numerical experiment is done on a parallel system with MIMD (Multiple Data, Multiple Instruction) structure and shared memory [22]. The language of parallel processing that is used in this paper is standard MPI [23].

Table 1 

The Specification of test bed.

The results of the procedure described above consist of different aspects. Its main purpose was to find the solution of a tridiagonal system of equations in a fast and efficient manner. The other achievements were computing the orthogonal and upper-triangular matrices which are useful computations in different fields. This entire procedure in a sequential form requires 𝑂(𝑛²) number of operations, where𝑛 is the number of equations. The complexity of the parallel algorithm is9𝐻²+26𝐻+32𝑛+130, (5.1) where 𝐻=𝑛/𝑝 is the ratio of the number of equations to the number of processors. Furthermore, as mentioned in Section 3, the number of zero elements raises as the number of blocks increases. This feature can be thought of as being on route to linear complexity.

Choosing a proper method to solve tridiagonal system of equations depends on different criteria, such as the order of complexity, accuracy and variety of achievements. For instance, the order of complexity of Gaussian elimination without pivoting is𝑂(𝑛), but it is not well suited for parallel processing [24]. On the other hand, in spite of favorable order of complexity of parallel iterative algorithms, their main focus is just on solution of the system with good enough accuracy but without any further outcomes.

The complexities of different 𝑄𝑅 factorization methods in sequential form are brought here as follows [25]:(1)Gram-Schmidt: (2/3)𝑛³+𝑂(𝑛²),(2)Householder: (4/3)𝑛³+𝑂(𝑛²),(3)Given: (8/3)𝑛³+𝑂(𝑛²),(4)Fast Given: (4/3)𝑛³+𝑂(𝑛²).

There are some important points to justify our algorithm. The first is the accuracy of obtained results. Obviously, the accuracy of direct solvers, that is, mostly exact, is better than that of the iterative ones. The second item is the variety of results, such as computing the orthogonal and upper-triangular matrices. The 𝑄𝑅 decomposition is usually used to solve the linear least squares problem and find the eigenvalues of system. The third is the scalability of the algorithm, which is a crucial parameter for success of a parallel algorithm, that is, the proposed algorithm can employ more processors to solve larger problems efficiently. Consequently, development of this efficient algorithm based on the 𝑄𝑅 factorization is valuable and practical. Some sequential and parallel 𝑄𝑅 factorization methods are investigated by Demmel et al. [26].

The obtained results are mostly exact. The accuracy of the algorithm was investigated using two methods. The first method was a norm of ‖𝑋‖ to evaluate error over all computed unknowns, and the second one was a norm of ‖𝑄^𝑇𝑄−𝐼‖ to determine the accuracy of the orthogonal matrix. The level of error for systems of order 𝑛 = 10 000, 20 000, 30 000, and 40 000 is shown in Table 2. The order of error in these cases is acceptable. The final results, that is, the computed unknown, became more accurate in larger matrices. Moreover, the level of orthogonality is high in all cases.

Table 2 

A norm of ‖𝑋‖ to evaluate error over all computed unknowns and a norm of ‖𝑄^𝑇𝑄−𝐼‖ to determine the accuracy of the orthogonal matrix.

Furthermore, this algorithm includes a new and interesting parallel feature. The parallel algorithm can decrease the computational cost and the requested memory on a single processor. The MPI can be used on a single processor to divide it into different parts. In other words, the parallel algorithm based on the MPI library can be executed on different parts of a processor in the same way as it is executed on a network of parallel CPUs. The enormous number of equations causes the processor to return the “insufficient virtual memory” error. Increasing the number of blocks in the matrix results in increasing the amount of zero elements and decreasing the size of the requested memory (physical and virtual). In Figure 7, the size of the requested memory for a matrix of order 𝑛 = 20 000 and 50 000 is shown. For instance, computation of a matrix of order 𝑛 = 50 000 was not possible on a specified processor with a single block. However, the computation was done on this processor with 4 blocks using 4.75 GB of memory. When the number of blocks was 16, the used memory decreased to 1.48 GB. The consumed time to find the solution of a system of equations of order 𝑛 = 20 000 and 𝑛 = 50 000 is presented in Figure 8. It is shown that when the number of blocks increased, time consumption was reduced. The rate of reduction nearly vanished when the number of blocks reached 16.

Figure 7 

The size of requested memory in computation by a single CPU.

Figure 8 

Time of computation by a single CPU.

The key issue in the parallel processing of a single application is speedup, especially its dependence on the number of processors used. Speedup is defined as a factor by which the execution time for the application changes:𝑆_𝑝=𝑇_𝑠𝑇_𝑝, (5.2) where, 𝑝, 𝑇_𝑠, and 𝑇_𝑝 are the number of processor, the executing time of the sequential algorithm, and the executing time of the parallel algorithm with 𝑝 processors, respectively.

On the other hand, the efficiency of parallel processing is a performance metric which is defined as𝐸_𝑝=𝑆𝑝. (5.3) It is a value, typically between zero and one, estimating how well utilized the processors are in solving the problem, compared to how much effort is wasted in communication and synchronization.

To investigate the ability of the parallel algorithm, the CPU time (𝑇), the speedup (𝑆_𝑝), and the efficiency (𝐸_𝑝) were considered in different cases (Table 3). The mentioned orders are 20 000, 40 000, 60 000, 80 000, and 100 000. The first column shows the computational time in a sequential case. A critical problem in solving a large system of equations using one processor, as mentioned before (e.g., order more than 50 000 in this study), is insufficient memory. To conquer this obstacle, the virtual memory was increased to 10 GB. Although using virtual memory can remove the “insufficient memory” error, it causes the speed of the computation to decrease, leading to inefficiency. The memory problem is completely solved in the parallel algorithm because each processor just needs to have the data of its block. For instance, if the order of matrix is 100 000, then in two processors the order of the data is 50 000, and in three processors it will reach to about 33 000. Therefore, the size of a huge system of equations does not create any computational limitation on the parallel cases.

Table 3 

The computation time, speed up, and efficiency of parallel solution of the systems of equations with different order.

The speedup and efficiency curves for the computation of a tridiagonal system of equations from orders 20 000 to 100 000 using different numbers of processors are presented in Figures 9–11. The comparison between computation times of different sizes of systems of equations using up to 16 processors is shown in Figure 9. The rate of time reduction was greater in larger systems. For instance, the time was reduced from 267 milliseconds to 22 milliseconds, that is, 11 times less than sequential form.

Figure 9 

The computation time of system of equation in different orders using up to 16 processors.

Figure 10 

The speedup of parallel algorithm.

Figure 11 

The efficiency of parallel algorithm.

The amount of speedup in different problem sizes is shown in Figure 10. When the order of the system was 20 000 and the number of processors was greater than 12, the speedup decreased. The reason for this inefficiency is the effect of network delay on computation time. On the other hand, in the larger problem sizes, this problem does not exist.

It is encouraging to note that, with an increase in problem size, the efficiency of the parallel algorithm increases. In Figure 11, the efficiency of different executions is presented. It is obvious that the rate of reduction in the efficiency of the parallel algorithm decreased by raising the problem size. This characteristic makes this parallel algorithm a practical subroutine in numerical computations.

In addition, the speedup of proposed algorithm is compared with parallel Cycle Reduction algorithm (Figure 12). Cycle Reduction is a well-known algorithm for solving system of equations [27]. The speedup of both algorithms is shown in two orders of system of equations, that is, 𝑛 = 200 000 and 400 000. Obviously, the speedup of proposed algorithm is more than that of Cycle Reduction and also would be greater by increasing number of equations and processors.

(a)

(b)

Figure 12 

The speedup comparison between proposed algorithm and Cycle Reduction in two different orders of system of equations: (a) 𝑛 = 200 000 and (b) 𝑛 = 400 000.

6. Conclusion

In this paper, a novel parallel algorithm based on the Gram-Schmidt method was presented to solve a tridiagonal linear system of equations in an efficient way. Solving linear systems of equations is a famous problem in different scientific fields, especially systems in the tridiagonal form. Therefore, the efficiency and accuracy of the algorithm can play a vital role in the entire solution procedure.

In this parallel algorithm, the coefficient matrix is partitioned according to the number of processors. After that, some columns in each block are switched. Therefore, the corresponding orthogonal and upper-triangular matrices are computed by each processor with a minimum request from the other partitions' data. Finally, the system of equations is solved almost exactly.

The proposed algorithm has some remarkable advantages, which include the following:(1)the ability to find almost exact solutions to tridiagonal linear systems of equations of different sizes, as well as the corresponding orthogonal and upper-triangular matrices;(2)considerable reduction of computation cost proportional to the number of processors and blocks;(3)increasing the efficiency of the algorithm by raising the size of the system of equations;(4)the significant ability of a single processor using the MPI library; the requested memory decreases by partitioning the matrix, and hence the computational speed increases.

The orders of the systems of equations that are presented in this paper range from 10 000 to 100 000. The results show that the suggested algorithm is scalable and completely practical.