BestSF: A Sparse Meta-Format for Optimizing SpMV on GPU BestSF: A Sparse Meta-Format for Optimizing SpMV on GPU

COO uses three arrays: one to store the nonzero elements of the matrix ($data$), and the two others to store the coordinates ($row\_index$ and $col\_index$). It is the most intuitive storage scheme for sparse matrices (Figure 1). The memory space needed for this format is $ size_{value} \times nnz + 2 \times size_{index} \times nnz$.

CSR uses three arrays: $data$, $col\_index$, and $ptr$. The array $data$ stores the nonzero elements. The integer array $col\_index$ is for storing the column indices of the nonzero elements of the sparse matrix. The last array, $ptr$, is used to store row pointers to the offset of each row in data (Figure 1). Due to its memory space efficiency ($ {size_{value} \times nnz + size_{index} \times (nnz+m+1)}$), CSR is the most popular format for storing sparse matrices.

BCSR divides the input matrix into blocks of $r \times c$ elements, and stores each non-empty block similar to the CSR format. If $r = c = 1$, the BCSR format is equivalent to the CSR format. An example of $2 \times 2$ block size is shown in Figure 1. BCSR format stores only one column index and one row pointer per block. Thus, the memory requirement is $size_{value} \times (r \times c \times k) + size_{index} \times ((n/r + 1) + k)$, where $k$ is the number of blocks. BCSR format can store fewer row pointers and column indices than CSR, but at the possible cost of filling in explicit zeros (padding).

The best block size that leads to the best performance depends on the sparsity of the input matrix (problem size, block density, etc.) [14]. To find the best block size for each input matrix, we considered a similar approach to [52, 53] based on the fill ratio estimation algorithm proposed in [54] which accurately measures the ratio of the number of stored values (including explicit zeros) to the number of non-zeros by sampling only a small number of rows of the matrix.

ELL uses two arrays: $data$ and $col\_index$. The array $data$ stores the nonzero elements. The integer array $col\_index$ stores the column indices of each nonzero element (Figure 1). If the dimension of the sparse matrix is $ n \times m$, each array has a size of $ n \times max$ with $max$ the maximum number of nonzero elements in a single row of the sparse matrix; thus, the memory space needed for this format is $ {size_{value} \times n \times max + size_{index} \times n \times max}$. Each row of $data$ stores the nonzero elements of the corresponding row of the sparse matrix, and each row from $col\_index$ stores the column indices of nonzero elements from the corresponding row of the sparse matrix. Zeros are added to every row of $data$ and $col\_index$ with the number of nonzero elements less than $max$ (zero padding). A given matrix fails to convert to ELL if the memory space needed after the zero padding exceeds the total global memory of the used GPU.

DIA uses two arrays: $data$ to store the nonzero elements, and offsets to store the offset of each diagonal from the central main diagonal which corresponds to the offset 0. DIA format is designed mainly for the sparse matrices for which most of the nonzero values are stored along diagonals. If some weak diagonals have very few nonzero elements, the DIA format may waste a large amount of memory space. In the example of the matrix A (Figure 1), DIA introduces a significant amount of zero padding, but it can be very suitable for regular sparse matrices arising from some regular discretization methods [6]. If $dia$ is the number of diagonals of the input sparse matrix with at least one nonzero element, then the memory needed for the DIA format is $ {size_{value} \times min(n,m) \times dia + size_{index} \times dia}$. Similar to ELL, a given matrix fails to convert to DIA if the memory space needed after the zero padding exceeds the total global memory of the used GPU.

HYB partitions the matrix into two parts, a “dense” part to be stored using ELL and a “very sparse” part to be stored using COO. Given an input sparse matrix, CUSP implementation computes a histogram of row sizes to find a threshold value $k$, where other implementations just take, as a threshold, the value of nonzero elements per row [47]. ELL is used for all the nonzero elements of the columns on the left side of $k$, whereas COO is used for the rest of the nonzero elements on the right of $k$ (Figure 1). The memory space needed for this format is $ size_{value} \times (n \times k + X) + size_{index} \times (n \times k + 2X)$ where $X$ is the number of nonzero values handled by the COO format.

To prove the importance of this problem, we show in the following the impact of selecting the best performing format on the overall performance of the SpMV kernel on GPU. To do so, we measured the performance of the SpMV kernel under the six considered sparse formats on a dataset¹ of 1,000 real-world sparse matrices from the SuiteSparse matrix collection [20]. All these sparse matrices fit into the global memory of the two GPUs used in our work (Table 1).

The pie chart in Figure 3 displays the distribution of the best performing sparse format on the 1,000 sparse matrices of our dataset. We can observe that, on the two GPUs, CSR covers an important percentage of the dataset (34% on Maxwell and 45% on Pascal). On the Maxwell GPU, BCSR comes in the second place followed by ELL and HYB. On the Pascal GPU, ELL comes in the second place followed by BCSR and HYB. On both GPUs, COO and DIA cover only a small percentage (less than 7%) of our dataset. Table 2 shows the average performance loss, on our dataset, if one of the sparse formats COO, CSR, BCSR, ELL, DIA, or HYB is used instead of the best performing format for each input matrix of our dataset. For example, the use of the HYB format introduces an average performance loss of more than 26% on the Maxwell GPU and 20% on the Pascal GPU. DIA and ELL introduce an important performance loss on the two GPUs. Note that given the zero padding needed in the construction of the sparse formats DIA and ELL, some sparse matrices fail to convert to these two sparse representations introducing 100% of performance loss. We can also notice that CSR has the best overall performance with less than 21.3% of performance loss on both GPUs which qualifies CSR as the “winner-takes-all” sparse format. However, as shown in Figure 4, it introduces a significant performance loss of more than 20% on an important number of sparse matrices of our dataset. All the previous observations clearly indicate that selecting the best format for each input matrix is very important for optimizing the overall performance of the SpMV kernel on GPU.

Fig. 3. Distribution of the best performing sparse format for the SpMV kernel on a dataset of 1,000 sparse matrices using two NVIDIA GPUs.

Fig. 4. Distribution of the matrices for which CSR is not the best performing sparse format on the intervals of performance loss of using CSR instead of the best performing sparse format for each matrix.

The metric used in our study to characterize the energy efficiency is performance per watt (FLOPS/W) or the number of floating point operations per joule (FLOP/Joule) as $ FLOPS/W = \frac{FLOP/Second}{Joules/Second} = FLOP/Joule$ [51]. For the SpMV kernel, the number of floating point operations is proportional to the number of nonzero elements of the sparse matrix, but the total energy (Joules) must be first measured in order to compute the performance per watt after that. We utilize the onboard GPU power sensors to obtain the instant power consumption of our kernels using the NVIDIA Management Library (NVML) interface which returns the power readings in milliwatts [43]. We followed a similar methodology as in [10, 39] to compute the total GPU energy consumed by our kernels.

We measured the energy efficiency of the SpMV kernel under the considered six sparse formats on our two GPUs for all the sparse matrices of our dataset. After that, we found the best sparse format in terms of energy efficiency for each one. Similar to Table 2, Table 3 shows the average energy efficiency loss, on our dataset, if one of the sparse formats COO, CSR, BCSR, ELL, DIA, or HYB is used instead of the most energy efficient format for each input matrix of our dataset. It is clear from Table 3 that there is no best for all sparse format in terms of energy efficiency. For example, using the CSR sparse format for all the matrices will introduce more than 22% of energy efficiency loss. Also, we noticed that, for 10% of the matrices of our dataset in Maxwell and almost 24% in Pascal, the most energy efficient sparse format is different from the best performing one. This observation joins what has been noticed in [16] and confirms that changing the program implementation (in our case the SpMV implementation using different sparse formats) may affect the performance and energy efficiency in a different way. However, we can notice in Figure 5, which shows the distribution of the matrices of our dataset on intervals of energy efficiency loss if the best sparse format in terms of performance is used instead of the most energy efficient one, that we have a significant energy efficiency loss (of more than 10%) only for a relatively small number of cases. All the previous observations suggest that using the best performing sparse format (GFLOPS) will also ensure the best energy efficiency (FLOPS/W) for a large number of sparse matrices.

Fig. 5. Distribution of the matrices of our dataset on the intervals of energy efficiency loss introduced by using the best performing sparse format instead of the most energy efficient sparse format.

As shown in Figure 6, the offline training of the pairwise models includes the following steps:

1. Constitute a dataset of sparse matrices that will be used in the learning and the testing of the classification system.
   
2. Extract some relevant sparsity features of the matrices that can capture at best the performance patterns of the considered sparse formats on GPU.
   
3. For all the matrices of the dataset, compute these features and run the SpMV kernel under each sparse format on the GPU to determine the real class of each matrix.
   
4. Select the most important features to use for training each pairwise model ($f_{i}$, $f_{j}$) and construct the pairwise models using a suitable machine learning algorithm. When training each model ($f_{i}$, $f_{j}$), each training point is weighted according to the performance difference between the sparse formats $f_{i}$ and $f_{j}$.
   

Given an input sparse matrix in a default sparse format $f^{(d)}$, the following steps are performed to predict its best sparse format:

1. Compute the sparsity features.
   
2. Use these features as inputs to the trained pairwise models. The predicted best sparse format, $f^{(p)}$, is then the one that has been predicted most often by the pairwise models.
   
3. Convert the considered sparse matrix to the predicted best sparse format ($f^{(d)} \rightarrow f^{(p)}$).
   

In our work, we used the implementation of the SpMV kernel under COO, CSR, DIA, ELL, and HYB available in the open source CUDA-based CUSP library [18]. Concerning the BCSR format, cuSPARSE [45] library is used. The dataset used includes 1,000 sparse matrices from the SuiteSparse matrix collection [20]. In the following sections, we provide all the details about the sparsity features used in our work, the learning process, and the evaluation of BestSF.

We chose to use SVMs as a training algorithm in building BestSF. SVMs, in general, have been identified in [30] as a promising machine learning method for the search algorithm selection problem in comparison with other machine learning methods. Also, SVMs deliver a good accuracy when the decision boundaries between different classes are highly nonlinear [17, 49], which we intuitively expect to be the case in the context of sparse format selection. Basically known as a binary classifier, several methods exist to extend SVMs for multiclass classification [27]:

• One-against-all approach: in the learning phase, it constructs $p$ SVM models where $p$ is the number of classes. The $n$th model is trained with all of the examples in the $n$th class with positive labels, and all other examples with negative labels. To find the predicted class for a new example, the decision functions of the previous models are evaluated and the new example is affected to the class with the highest value of the decision function.
  
• One-against-one approach: in the learning phase, this method performs a pairwise classification by constructing $p(p - 1)/2$ classifiers where each one is trained on data from two classes. Given a new example, it is affected to the class that has been voted most often by the previously trained modes. It has been shown in [27] that the one-against-one approach is more suitable for practical use than the one-against-all. This approach was used in our previous work [8].
  

The main drawback of the previous methods based on standard SVMs is that the training points are all given equal importance in the training phase and does not allow relative importance of training points to be considered. To solve this problem, we used the binary WSVM [60] to build $p \times (p-1) / 2$ WSVM pairwise models ($p$ is the number of sparse formats). For training the model $(f_{i}, f_{j})$, we attribute a weight, $w_{s}$, to each training point ($s$) according to the performance difference between the corresponding pair of sparse formats as shown in Equation (1). $P(A_{s},f_{i})$ and $ P(A_{s},f_{j})$ in Equation (1) represent the performance of the SpMV kernel using the matrix $A_{s}$ under the sparse formats $f_{i}$, $f_{j}$ respectively. $w_{s}(f_{i},f_{j})$ depends linearly on the performance improvement between the two sparse formats $f_{i}$ and $f_{j}$. Given a new example (sparse matrix in our case), we used a similar voting strategy to the one used in the standard one-against-one strategy: if the model $(f_{i}, f_{j})$ selects $f_{i}$, then the score of $f_{i}$ is increased by one. Otherwise, the score of $f_{j}$ is increased by one. We repeat this process with all $p \times (p-1) / 2$ pairwise models. Then, the new example is affected to the class with the largest score.

We used the implementation of WSVM available in the open source LibSVM library [13]. We kept using the standard C-SVM algorithm and selected the Gaussian Radial Basis Function (RBF) kernel. We used the cross-validation technique as suggested in [26] for selecting the two parameters $ C$ and $ gamma$. The parameter $ \alpha$ in Equation (1) is set to $ \alpha = 100.0$ (we achieved our best results with this value). A very important step is scaling the data points before the training/testing process. LibSVM comes with a default linear scaling scheme that linearly scales each attribute to the range $[-1,1]$ or $ [0,1]$. However, we found that using a logarithmic transformation ($ {Log_{10}(x)}$) of our data points instead of the linear scaling improves both the accuracy of classification and the performance of BestSF.

Our feature set includes 11 sparsity features of the sparse matrices as shown in Table 4. Besides the three-dimensional features $n$, $m$, and $nnz$, that characterize the volume of work needed to compute the SpMV kernel, not all of the other features are important for training each pairwise model. In other words, we have to select the most important features for each pairwise model. To do so, we used a wrapper-based approach [28] that runs the learning process with different feature subsets of our initial feature set (our 11 features) to determine the best subset to train each pairwise model. Given that our initial feature set includes only 11 features, we used the exhaustive search as a base for our wrapper feature selection, but other search strategies can be used if more features are included and the search space becomes very huge. The features subsets in Table 5 are obtained using our dataset on our GPUs. These features subsets may slightly change if we use other GPUs or other datasets. This is why BestSF includes the feature selection as part of the learning phase as shown in Figure 6.

First, we run the SpMV kernel under the considered sparse formats for all the 1,000 sparse matrices of our dataset, and each matrix is attributed to its class according to its best sparse format. Then, we randomly select 80% of our dataset (800 sparse matrices) to constitute the training set that will be used to train each of the pairwise models. The remaining 20% (200 sparse matrices) are used as a testing set. This 80-20 splitting of our dataset is repeated five times to constitute five different learning/testing sets. We run the training process to train all the pairwise models using WSVN as explained above. Then, the testing set is used to evaluate the accuracy, performance, and energy efficiency of BestSF.

Given a testing set of $N$ unseen sparse matrices ($N = 200$ for our case) and if $N_{+}$ is the number of matrices correctly classified, then the accuracy is calculated as $ 100 \times N_{+}/ N$. Table 6 shows that BestSF achieved from 83.4% to 85.3% of accuracy between the two used GPUs.

We calculate the average Performance Loss Under Best (PLUB). This metric measures how far the performance of BestSF is from the performance of an ideal classifier (with 100% of accuracy). If we consider a testing set with $ N$ sparse matrices, $ f^{(r)}_{k}$ is the real best performing sparse format for each sparse matrix $ A_{k}$ from the testing set, then PLUB is given in Equation (2). We also calculate PGO($f_{i}$), the average Performance Gain of using BestSF Over using the sparse format $ f_{i}$ for all the testing matrices as shown in Equation (3). Table 6 shows the PLUB values obtained on our two GPUs and also the values of PGO(COO), PGO(CSR), PGO(BCSR), and PGO(HYB). Note that we did not report the values of PGO(ELL) and PGO(DIA) because some testing matrices fail to convert to the ELL and DIA formats, which makes Equation (3) not applicable for all the testing points.

We summarize our observations in the following points:

• On both GPUs, the PLUB is less than 3%, which means that by using BestSF to select the best performing sparse format we reached more than 97% of the maximum performance possible. By checking the misclassified testing points, we found that most of the misclassifications occur when BestSF selects the wrong sparse format but with a very close performance to the optimal one. This explains why we have a near-optimal performance even if the accuracy of classification did not exceed 86% on both GPUs.
  
• Compared with just using one format (COO, CSR, BCSR, or HYB) for all the testing sparse matrices, BestSF registered an average performance improvement ranging from 41% to 94.3% confirming that using BestSF is better, in terms of average performance, than using just one sparse format for all the sparse matrices of our testing sets.
  

We recall that BestSF is trained to select the best performing (GFLOPS) sparse format and not the most energy efficient format. We used the metric ELUB to measure the average energy efficiency loss introduced by using BestSF under using the most energy efficient sparse format for all the inputs. If we consider $N$ testing matrices, $EE(A_{k},f_{i})$ is the GPU energy efficiency of the SpMV kernel using the sparse format $f_{i}$, $f_{k}^{(e)}$ is the most energy efficient sparse format for each sparse matrix $ A_{k}$ of the testing set, then ELUB is given in Equation (4). We recall that $f_{k}^{(p)}$ is the predicted sparse format for the matrix $ A_{k}$. We also calculate EGO($f_{i}$), the average Energy efficiency Gain of using BestSF Over using the sparse format $ f_{i}$ for all the testing matrices as shown in Equation (5). Table 7 shows the ELUB values obtained on our two GPUs and also the values of EGO(COO), EGO(CSR), EGO(BCSR), and EGO(HYB).

We can make the following observations:

• On the two used GPUs, the ELUB is less than 4% meaning that using BestSF achieved more than 96% of the optimal energy efficiency possible.
  
• Compared with just using one sparse format (COO, CSR, BCSR, or HYB) for all the testing sparse matrices, BestSF registered an average energy efficiency improvement ranging from 46.1% to 96.2% confirming that using BestSF is better, in terms of average energy efficiency, than using just one sparse format for all the input matrices.
  

In the following, we first compare the pairwise weighted classification approach used in BestSF to two other learning approaches, standard multiclass classification and decision trees, previously used for solving the problem of sparse format selection [8, 40, 41, 47]. After that, we use a set of commonly evaluated large sparse matrices to compare the performance of BestSF with the performance of CSR5 [35] and merge-based CSR [36].

In [47], decision trees were used for training a classification system for selecting the best sparse format from CSR, ELL, and HYB. We showed in our previous work [8], that using multi-class SVMs with very simple to compute sparsity features allows better accuracy and selection quality. In [40, 41], the authors proposed Nitro to select the best variant of different kernels (CG solvers, BFS, Histogram, Sort, and SpMV) based on SVMs. Concerning the SpMV kernel, Nitro considers only the three sparse formats CSR, ELL, and DIA. To improve the performance of the SpMV kernel, BestSF includes six sparse formats COO, CSR, BCSR, ELL, DIA, and HYB. In addition, to improve the quality of selection, it relies on a pairwise weighted classification approach based on WSVM as explained previously. Another important aspect that makes the difference between BestSF and the previous solutions is the set of used sparsity features. Table 8 exposes this aspect and also shows the accuracy of classification, PLUB, and ELUB obtained with using the standard multi-class classification and decision trees (BFtree). Comparing to Table 6 and Table 7, we can notice that the weighted pairwise classification, used in BestSF, outperforms both the standard one-vs-one multi-class SVMs and BFTrees in terms of accuracy, average PLUB, and average ELUB values on the two used GPUs. Also, using our sparsity features improves all the previous metrics. In addition, Figure 7 shows the distribution of the misclassifications on intervals of PLUB values obtained with different learning approaches. We can see that using the weighted pairwise classification reduces the number of misclassifications with more than 10% of performance loss (PLUB > 10%). All the previous observations confirm that the weighted pairwise classification based on WSVM, used in BestSF, is more suitable for the sparse format selection problem.

Fig. 7. Distribution of misclassifications on intervals of PLUB.

Figure 8 compares the performance and energy efficiency of the SpMV kernel under BestSF, across commonly used sparse matrices in evaluating different sparse formats [6, 35], to the performance obtained with using merge-based CSR [36] and CSR5 [35]. Unlike many sparse formats that perform very good on some sparse matrices and very bad on others, both CSR5 and merge-based CSR have been designed to deliver consistently high performance on a wide range of sparsity patterns (structured and unstructured). The merge-based CSR is a new implementation of the SpMV kernel using the CSR format based on a fine-grained merge-based parallel decomposition to have a well-balanced execution. CSR5 is a new variation of the CSR format and the CSR5-based SpMV kernel is implemented based on a new low-overhead segmented sum algorithm. It showed a good performance in comparaison to HYB, BRC (Blocked Row-Column [5]), and ACSR (Adaptive-CSR [4]) on both structured and unstructured matrices. We can observe in Figure 8 that BestSF, by selecting the right sparse format (from COO, CSR, ELL, DIA, HYB, and BCSR) for each sparse matrix, achieved an important performance improvement for most of the considered sparse matrices. We conclude that selecting between different existing formats is a very promising approach to have a consistently high performance across different sparsity patterns.

Fig. 8. Comparing BestSF to CSR5 [35] and merge-based CSR [36] on commonly evaluated large sparse matrices.

As shown in Figure 6, BestSF includes an online decision-making phase in which given an input sparse matrix in its default sparse format $f^{(d)}$, the sparsity features are first calculated, the pairwise models are used to predict the best sparse format $f^{(p)}$, and finally the input sparse matrix is converted to the sparse format $f^{(p)}$ ($f^{(d)} \rightarrow f^{(p)}$). Thus, the practical effectiveness of BestSF depends on the time overhead to perform this online decision-making phase. Note that all the steps of the online decision-making phase are executed on the CPU side (Intel i7-3770 running at 3.4GHz for our machine) using sequential implementations (we implemented our algorithms for computing the sparsity features, and used CUSP and cuSPARSE routines for the sparse format conversion).

Table 9 shows the sparse matrices used to evaluate the time overhead of BestSF. The sparse matrices of G1 and G2 are unstructured sparse matrices selected from the SuiteSparse matrix collection [20]. G1 matrices have been suggested in [32] and most of G2 matrices have been used in [11, 12]. The two sparse matrices of G3 are structured matrices resulting from a 3D finite difference discretization using the standard seven-point stencil. The matrices G2+G3 have the particularity of being square, symmetric, and positive-definite, and they are also used to evaluate the preconditioned conjugate gradient iterative solver based on BestSF in Section 8.

If the time needed for the prediction, including both features extraction and the runtime of the models, is $ p\_time$ and the time needed for sparse format conversion is $ c\_time$, then the total time overhead is $ tot\_time = p\_time + c\_time$. If the predicted sparse format $f^{(p)}$ is the same as the default sparse format $f^{(d)}$ (CSR in our case), then there is no format conversion needed and $ c\_time$ equals to zero in this case. Table 10 shows both $ p\_time$ and $ c\_time$ and also $tot\_time$ for all the sparse matrices of Table 9 for our two GPUs. As suggested in [32], the $tot\_time$ is also presented in terms of CSR-SpMV time on GPU. For example, an overhead of 30 $\times$ CSR-SpMV means that the time overhead for the online decision-making phase takes about 30 times as long as the actual GPU SpMV time using the standard CSR representation. We refer to the CSR format as it is the most used format in different application domains. As shown in Table 10, the overhead in terms of CSR-SpMV, is ranging from 21 to 152 on the Maxwell GPU and from 17 to 144 on the Pascal GPU. We can notice that, except the matrices for which $ c\_time$ is zero ($f^{(p)}$ = CSR), the sparse format conversion is the main source of overhead. Note that the need for format conversion is not specific to only BestSF. For example, if we use ELL or HYB for all the input matrices (without using BestSF), putting a sparse matrix in these sparse formats usually needs to be done from the basic COO or CSR formats [21]. Note that $p\_time$ is largely dominated by the time needed for features extraction. Using the trained WSVM models is relatively fast because it mainly consists of simple dot products between the features vector (representing the input matrix) and the support vectors of the models. The measured runtime of the models is $ 0.048$ms for the Maxwell GPU and $ 0.056$ms for the Pascal GPU, which represents less than 1% of the total overhead for all the testing matrices (Table 10).

Given the overhead of the online decision-making phase, we conclude that BestSF is more suitable for applications in which the SpMV kernel is computed many times with the same sparse matrix. This is the case of a wide range of applications using iterative solvers for solving large linear systems ($Ax = b$) and eigenvalue problems ($Ax = \lambda x$) where the SpMV kernel is usually executed hundreds if not thousands of times with the same matrix $A$ [6, 31].

Conjugate Gradient is one of the most popular iterative methods for solving large systems of linear equations $ A x = b$ where $x$ is an unknown vector, $b$ is a known vector, and $A$ is a known, square, symmetric, positive-definite matrix [46, 48]. This method is generally applicable to large sparse systems for which direct methods are not efficient. This kind of sparse system arises very often when numerically solving partial differential equations and optimization problems. For unstructured sparse matrices, generally characterized by a large condition number (the ratio of the largest to smallest eigenvalue), the convergence of the CG method may need a large number of iterations. Preconditioning is usually used as a technique for improving the condition number of the input matrix and reducing the number of iterations required for convergence. However, an efficient preconditioner is not only the preconditioner that considerably reduces the number of iterations needed for the convergence but also the preconditioner for which the preconditioning matrix is easy to generate, and the preconditioning operation is easy to apply [33, 46]. Different preconditioners can be used with the PCG formulation shown in Figure 9. In this work, we consider the Jacobi (diagonal) preconditioner because of its embarrassingly parallel structure [22]. Other preconditioners, like the Symmetric Successive Over Relaxation (SSOR) and the Incomplete LU (ILU), offer better convergence rates, but applying these preconditioners comes at the price of expensive sparse triangular equations that can become a bottleneck on parallel architectures [1, 2, 34].

Fig. 9. Preconditioned Conjugate Gradient (PCG) algorithm.

In this work, we use the formulation of the PCG algorithm shown in Figure 9. In each iteration (lines 8–17), the algorithm uses two dot product (DOT) operations (lines 9 and 14), two linear combinations of vectors (AXPY) in lines 10 and 11 and one AYPX operation in line 16 which consist of two successive operations: scaling a vector by a scalar (SCAL) and linear combination of vectors (AXPY). As we are considering the diagonal preconditioner, the preconditioning operation (line 12) consists of an element-wise multiplication (XMY) of the diagonal $d$ of the matrix A ($ M = d^{-1}$) by the vector $r$. The SpMV operation (line 8) is by far the most expensive operation among all the operations used by the Jacobi-based PCG algorithm.

In our implementation, the NVIDIA's cuBLAS library [44] was used for the kernels DOT, AXPY, and AYPX (by combining SCAL and AXPY). For the XMY kernel, the thrust library [7] was used. Finally, for the SpMV kernel, we used the implementation of this kernel under different sparse formats COO, CSR, DIA, ELL, and HYB available in CUSP library and BCSR available in cuSPARSE.

To evaluate the gain of using BestSF for the SpMV kernel inside the PCG solver, we measured the total execution time (including the CPU-GPU data transfer) of the PCG algorithm with different sparse formats COO, CSR, BCSR, ELL, DIA, and HYB. We also used CSR5 and merge-based CSR to be able to compare BestSF with exterior sparse formats. We used 14 square, symmetric, positive-definite sparse matrices (Table 9, G2 and G3). We recall that the two matrices poisson192 and poisson256 (G3) are resulting from a 3D finite difference discretization using the standard seven-point stencil. They are banded and diagonally dominant structured with a relatively low condition number. We included them in the testing matrices to demonstrate the case of having a relatively small number of iterations. Note that all 14 testing matrices have not been included in the learning dataset used for training BestSF. In our experiments, the tolerance factor $ \epsilon$ is set to $ \epsilon = 10^{-6}$ and double precision is used.

Table 11 shows the execution time needed for the convergence for PCG-BestSF with our 14 sparse matrices and also the performance gain (PGO($f_{i}$)) of using BestSF over using each of its constituent sparse formats (COO, CSR, BCSR, DIA, ELL, and HYB) and also over using merge-based CSR [36] and CSR5 [35]. Negative values of PGO($f_{i}$) indicate a performance loss under using the sparse format $f_{i}$. Note that the execution time reported in Table 11 represents the total execution time including the CPU-GPU data transfer, sparse format conversion, and format prediction time. For all 14 testing sparse matrices, we registered only one misclassification on each GPU (G3_circuit with 1.5% and 1.4% of performance loss on Maxwell and Pascal, respectively). BestSF was able to make the best choice for the matrices where there is a relatively important performance difference between different sparse formats like in crankseg_2 and Bump_2911. Other than the performance loss introduced by a misclassification, there is a performance loss under the sparse format that has been selected for each matrix due to the time overhead discussed previously in Section 7.5. This performance loss does not exceed 3% on both GPUs even for the matrices with a relatively small number of iterations like poisson192 and poisson256. An important performance improvement is also registered in comparison with CSR5 and merge-based CSR for most of the considered matrices.

Table 12 shows the GPU total energy consumption for PCG-BestSF with our 14 testing sparse matrices. Note that the energy consumption reported in Table 12 does not include the energy consumed by the CPU (needed for controlling the GPU). We also reported the energy gain (EGO($f_{i}$)) of using BestSF over using its constituent sparse formats (COO, CSR, BCSR, DIA, ELL, and HYB) and also CSR5 and merge-based CSR. Negative values of EGO($f_{i}$) indicate an energy loss under using the sparse format $f_{i}$. EGO($f_{i}$) is reported 0 when the sparse format selected by BestSF ($ f^{(p)}$) is also the most energy efficient sparse format ($f^{(e)}$). We observe that BestSF achieved an important energy efficiency gain over using the popular CSR format (EGO(CSR)) for matrices like F1_639 and Bump_2911. We can make similar observations concerning the other sparse formats. On both GPUs, we have only one matrix for which the predicted sparse format by BestSF is different from the most energy efficient format (G3_circuit with 10.6% and 15.7% of energy efficiency loss, respectively).

All the previous observations prove the practical effectiveness of BestSF, especially for unstructured sparse matrices for which the winner-takes-all approach to the format selection problem leads to an important overall performance and energy efficiency loss.