Performance Analysis of Effective Symbolic Methods for Solving Band Matrix SLAEs

This paper presents an experimental performance study of implementations of three symbolic algorithms for solving band matrix systems of linear algebraic equations with heptadiagonal, pentadiagonal, and tridiagonal coefficient matrices. The only assumption on the coefficient matrix in order for the algorithms to be stable is nonsingularity. These algorithms are implemented using the GiNaC library of C++ and the SymPy library of Python, considering five different data storing classes. Performance analysis of the implementations is done using the high-performance computing (HPC) platforms"HybriLIT"and"Avitohol". The experimental setup and the results from the conducted computations on the individual computer systems are presented and discussed. An analysis of the three algorithms is performed.


Introduction
Systems of linear algebraic equations (SLAEs) with heptadiagonal (HD), pentadiagonal (PD) and tridiagonal (TD) coefficient matrices may arise after many different scientific and engineering problems, as well as problems of the computational linear algebra where finding the solution of a SLAE is considered to be one of the most important problems. On the other hand, special matrix's characteristics like diagonal dominance, positive definiteness, etc. are not always feasible. The latter two points explain why there is a need of methods for solving of SLAEs which take into account the band structure of the matrices and do not have any other special requirements to them. One possible approach to this problem is the symbolic algorithms. An overview of some of the symbolic algorithms which exist in the literature is done by us in [1]. What is common for all of them, is that they are implemented using Computer Algebra Systems (CASs) such as Maple, Mathematica, and Matlab.
Three direct symbolic algorithms for solving a SLAE based on LU decomposition are considered: for HD matrices (see [2] and [1]) -SHDM, for PD [3] -SPDM, and TD matrices [4] -STDM. The only assumption on the coefficient matrix is nonsingularity.
The choice of algorithms for solving problems of the computational linear algebra is crucial for the programs' effectiveness, especially when these problems are with a big dimension. In that case they require the use of supercomputers and computer clusters for the solution to be obtained in a reasonable amount of time. Another important choice that has to be made and that influences the programs' performance is what programming language to be used for the algorithms' implementations. Here, we are going to focus on two of the most popular programming languages for scientific computations, namely C++ and Python. The aim of this paper, which is a logical continuation of [5] and [6], is to investigate the performance characteristics of the considered serial methods with their two implementations being executed on modern computer clusters.

Computational Experiments
Computations were held on the basis of the heterogeneous computational platform "Hy-briLIT" (1142 TFlops/s for single precision and 550 TFlops/s for double precision) [7] at the Laboratory of Information Technologies of the Joint Institute for Nuclear Research in the town of science Dubna, Russia, and on the cluster computer system "Avitohol" (412.3 TFlops/s for double precision) [8] at the Advanced Computing and Data Centre of the Institute of Information and Communication Technologies of the Bulgarian Academy of Sciences in Sofia, Bulgaria. The latter has been ranked among the TOP500 list (https: //www.top500.org) twice -being 332nd in June 2015 and 388th in November 2015. Table 1 sums up some basic information about hardware on the two computer systems, including models of processors, processors' base frequency, and amount of cache memory (Smart-Cache). For even more information, visit: https://ark.intel.com/compare/75281,75269. Table 2 summarizes the basic information about the compilers and libraries used on the two computer systems. The three algorithms are implemented using the GiNaC library [9] of C++ [10] (the projects are built with the help of CMake [11]), and SymPy [12] library of Python [13] (using Anaconda distribution [14]). The reason why optimization -O0 was used is that the GiNaC library is already optimized and any further attempts to optimize the code give worse results. Although we use Py 2.7, one should note that the implementations are fully compatible with Py 3.6 as well. It is worth mentioning that optimization attempts like autowrap and numba do not give better execution times. The former because of overhead, the second one because it cannot optimize further than what is already done.

Experimental Results and Analysis of the Algorithms
During our experiments wall-clock times were collected and the average time from multiple runs is reported. For that purpose, we use the function now() of the class std::chrono::high_resolution_clock for C++ (requires at least standard c++11), and timeit.default_timer() for Python. Both of them provide the best clock rate available on the platform. Five different classes for data storing are testedlst and matrix of GiNaC; Matrix (variable-size), Matrix (fixed-size), and MutableDenseNDimArray of SymPy. The first and the third are intended for variable-size storing, while the others -for fixed-size one (the SymPy's Matrix could be both). Further we are going to denote the implementations of the three algorithms, using these 5 data storing classes as Impl. i, i = 1, . . . 5. The notation is as follows: SXDM stands for symbolic X method, X = HD, PD, TD. The achieved computational times from solving a SLAE are summarized in Tables 3-7.  Table 4: Results from solving a SLAE on the two clusters applying the second implementation.
GiNaC ("matrix" -fixed-size) := Impl.    Using the following formula (k -unknown coefficient of proportionality, t i -time, N ithe matrix's number of rows, i = 1, 2): the order of growth of execution time α for all the five implementations was estimated (see Table 8).
Remark: The number of needed operations for Gaussian elimination so as a HD SLAE to be transformed into a PD one is 35N − 122, while the number of needed operations so as a PD SLAE to be transformed into a TD one is 23N − 52, where N is the matrix's number of rows. This observation is relevant to Figure 1

Discussion and Conclusions
Direct comparison between the five implementations shows that the best GiNaC one is Impl. 2, while the best SymPy one is Impl. 5 (see Figure 2). Expectedly, the GiNaC implementations of the three algorithms yield a much better computation time in comparison with the respected SymPy implementations with the difference being between one and two orders of magnitude. It must be mentioned that the matrix class in SymPy is a subclass of the ndarray. This means that every call on a matrix object requires a few extra Python calls. This usually leads to a little slower performance although the difference is negligible. Figures 3 and 4 depict a comparison of the execution time of the best GiNaC and SymPy implementations on the two clusters. As one can see, "HybriLIT" behaves better than "Avitohol" with the difference being bigger when the SymPy is of interest. It is obvious that the implementations which rely on variable-size storing classes (that are Impl. 2 and Impl.3) were found to be much slower than the ones which use fixed-size storing classes, but while this was expected, it does not belittle their importance since there is a class of problems where the size of the matrix is unknown at runtime.
As it can be seen in Table 9, three of the implementations have a linear time growth, while the other two have a quadratic trend. Theoretically, all the algorithms introduced here have   a linear complexity. This means that even if theoretically the execution time has to grow as O(N), the two implementations which rely on variable-size storing classes show time growth as O(N 2 ), with Impl. 1 even having α > 2. Hence, the variable-size of the storing class changes the order of complexity with one order of magnitude. The choice of a programming language depends on a lot of factors. However, in the context of symbolic computations for solving a SLAE with a band coefficient matrix (length of band equal to 3, 5 or 7) among the options suggested in this work, we can note the following: Python is easier to learn and easier and faster to prototype (no need for memory management, build systems, compilers, etc.), on the other hand, C++ has better performance and occupies less memory, but requires much more attention to bookkeeping and storage details.