\{\texttt{lilia.ziane\_khoja},~\texttt{raphael.couturier},~\texttt{arnaud.giersch},~\texttt{jacques.bahi}\}\texttt{@univ-fcomte.fr}
}
+\newcommand{\Iter}{\mathit{iter}}
+\newcommand{\Max}{\mathit{max}}
+\newcommand{\Offset}{\mathit{offset}}
+\newcommand{\Prec}{\mathit{prec}}
+\newcommand{\Ratio}{\mathit{Ratio}}
+\newcommand{\Time}[1]{\mathit{Time}_\mathit{#1}}
+\newcommand{\Wavg}{W_{\mathit{avg}}}
+
\begin{document}
\maketitle
\begin{abstract}
In this paper, we aim at exploiting the power computing of a GPU cluster for solving large sparse
linear systems. We implement the parallel algorithm of the GMRES iterative method using the CUDA
-programming language and the MPI parallel environment. The experiments shows that a GPU cluster
+programming language and the MPI parallel environment. The experiments show that a GPU cluster
is more efficient than a CPU cluster. In order to optimize the performances, we use a compressed
storage format for the sparse vectors and the hypergraph partitioning. These solutions improve
the spatial and temporal localization of the shared data between the computing nodes of the GPU
%%% END %%%
\begin{algorithm}[!h]
+ \newcommand{\Convergence}{\mathit{convergence}}
+ \newcommand{\False}{\mathit{false}}
+ \newcommand{\True}{\mathit{true}}
%\SetAlgoLined
- \Entree{$A$ (matrix), $b$ (vector), $M$ (preconditioning matrix),
-$x_{0}$ (initial guess), $\varepsilon$ (tolerance threshold), $max$ (maximum number of iterations),
+ \KwIn{$A$ (matrix), $b$ (vector), $M$ (preconditioning matrix),
+$x_{0}$ (initial guess), $\varepsilon$ (tolerance threshold), $\Max$ (maximum number of iterations),
$m$ (number of iterations of the Arnoldi process)}
- \Sortie{$x$ (solution vector)}
+ \KwOut{$x$ (solution vector)}
\BlankLine
$r_{0} \leftarrow M^{-1}(b - Ax_{0})$\;
$\beta \leftarrow \|r_{0}\|_{2}$\;
$\alpha \leftarrow \|M^{-1}b\|_{2}$\;
- $convergence \leftarrow false$\;
+ $\Convergence \leftarrow \False$\;
$k \leftarrow 1$\;
\BlankLine
- \While{$(\neg convergence)$}{
+ \While{$(\neg \Convergence)$}{
$v_{1} \leftarrow r_{0} / \beta$\;
\For{$j=1$ {\bf to} $m$}{
$w_{j} \leftarrow M^{-1}Av_{j}$\;
}
\BlankLine
Put $V_{m}=\{v_{j}\}_{1\leq j \leq m}$ and $\bar{H}_{m}=(h_{i,j})$ Hessenberg matrix of order $(m+1)\times m$\;
- Solve the least-squares problem of size $m$: $\underset{y\in\mathbb{R}^{m}}{min}\|\beta e_{1}-\bar{H}_{m}y\|_{2}$\;
+ Solve the least-squares problem of size $m$: $\displaystyle\min_{y\in\mathbb{R}^{m}} \|\beta e_{1}-\bar{H}_{m}y\|_{2}$\;
\BlankLine
$x_{m} \leftarrow x_{0} + V_{m}y$\;
$r_{m} \leftarrow M^{-1}(b-Ax_{m})$\;
$\beta \leftarrow \|r_{m}\|_{2}$\;
\BlankLine
- \eIf{$(\frac{\beta}{\alpha}<\varepsilon)$ {\bf or} $(k\geq max)$}{
- $convergence \leftarrow true$\;
+ \eIf{$(\frac{\beta}{\alpha}<\varepsilon)$ {\bf or} $(k\geq \Max)$}{
+ $\Convergence \leftarrow \True$\;
}{
$x_{0} \leftarrow x_{m}$\;
$r_{0} \leftarrow r_{m}$\;
It is quite expensive for large size matrices in terms of execution time and memory space. In
addition, it performs irregular memory accesses to read the nonzero values of the sparse matrix,
implying non coalescent accesses to the GPU global memory which slow down the performances of
-the GPUs. So we use the HYB kernel developed and optimized by Nvidia~\cite{CUSP} which gives on
+the GPUs. So we use the HYB kernel developed and optimized by NVIDIA~\cite{CUSP} which gives on
average the best performance in sparse matrix-vector multiplications on GPUs~\cite{Bel09}. The
HYB (Hybrid) storage format is the combination of two sparse storage formats: Ellpack format
(ELL) and Coordinate format (COO). It stores a typical number of nonzero values per row in ELL
\label{fig:02}
\end{figure}
-Linux cluster version 2.6.18 OS is installed on the six machines. The C programming language is used for
+Scientific Linux 5.10, with Linux version 2.6.18, is installed on the six machines. The C programming language is used for
coding the GMRES algorithm on both the CPU and the GPU versions. CUDA version 4.0~\cite{ref19} is used for programming
the GPUs, using CUBLAS library~\cite{ref37} to deal with the functions operating on vectors. Finally, MPI routines
of OpenMPI 1.3.3 are used to carry out the communication between the CPU cores.
The experiments are done on linear systems associated to sparse matrices chosen from the Davis collection of the
-university of Florida~\cite{Dav97}. They are matrices arising in real-world applications. Table~\ref{tab:01} shows
+University of Florida~\cite{Dav97}. They are matrices arising in real-world applications. Table~\ref{tab:01} shows
the main characteristics of these sparse matrices and Figure~\ref{fig:03} shows their sparse structures. For
each matrix, we give the number of rows (column~$3$ in Table~\ref{tab:01}), the number of nonzero values (column~$4$)
and the bandwidth (column~$5$).
All the experiments are performed on double-precision data. The parameters of the parallel
GMRES algorithm are as follows: the tolerance threshold $\varepsilon=10^{-12}$, the maximum
-number of iterations $max=500$, the Arnoldi process is limited to $m=16$ iterations, the elements
+number of iterations $\Max=500$, the Arnoldi process is limited to $m=16$ iterations, the elements
of the guess solution $x_0$ is initialized to $0$ and those of the right-hand side vector are
-initialized to $1$. For simplicity sake, we chose the matrix preconditioning $M$ as the
+initialized to $1$. For simplicity's sake, we chose the matrix preconditioning $M$ as the
main diagonal of the sparse matrix $A$. Indeed, it allows us to easily compute the required inverse
matrix $M^{-1}$ and it provides relatively good preconditioning in most cases. Finally, we set
the size of a thread-block in GPUs to $512$ threads.
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
-Matrix & $Time_{cpu}$ & $Time_{gpu}$ & $\tau$ & \# $iter$ & $prec$ & $\Delta$ \\ \hline \hline
+Matrix & $\Time{cpu}$ & $\Time{gpu}$ & $\tau$ & \# $\Iter$ & $\Prec$ & $\Delta$ \\ \hline \hline
2cubes\_sphere & 0.234 s & 0.124 s & 1.88 & 21 & 2.10e-14 & 3.47e-18 \\
ecology2 & 0.076 s & 0.035 s & 2.15 & 21 & 4.30e-13 & 4.38e-15 \\
finan512 & 0.073 s & 0.052 s & 1.40 & 17 & 3.21e-12 & 5.00e-16 \\
In Table~\ref{tab:02}, we give the performances of the parallel GMRES algorithm for solving the linear
systems associated to the sparse matrices shown in Table~\ref{tab:01}. The second and third columns show
the execution times in seconds obtained on a cluster of 24 CPU cores and on a cluster of 12 GPUs, respectively.
-The fourth column shows the ratio $\tau$ between the CPU time $Time_{cpu}$ and the GPU time $Time_{gpu}$ that
+The fourth column shows the ratio $\tau$ between the CPU time $\Time{cpu}$ and the GPU time $\Time{gpu}$ that
is computed as follows:
\begin{equation}
- \tau = \frac{Time_{cpu}}{Time_{gpu}}.
+ \tau = \frac{\Time{cpu}}{\Time{gpu}}.
\end{equation}
From these ratios, we can notice that the use of many GPUs is not interesting to solve small sparse linear
systems. Solving these sparse linear systems on a cluster of 12 GPUs is as fast as on a cluster of 24 CPU
cores. Indeed, the small sizes of the sparse matrices do not allow to maximize the utilization of the GPU
cores of the cluster. The fifth, sixth and seventh columns show, respectively, the number of iterations performed
-by the parallel GMRES algorithm to converge, the precision of the solution, $prec$, computed on the GPU
+by the parallel GMRES algorithm to converge, the precision of the solution, $\Prec$, computed on the GPU
cluster and the difference, $\Delta$, between the solutions computed on the GPU and the GPU clusters. The
last two parameters are used to validate the results obtained on the GPU cluster and they are computed as
follows:
\begin{equation}
-\begin{array}{c}
- prec = \|M^{-1}(b-Ax^{gpu})\|_{\infty}, \\
- \Delta = \|x^{cpu}-x^{gpu}\|_{\infty},
-\end{array}
+ \begin{aligned}
+ \Prec &= \|M^{-1}(b-Ax^{gpu})\|_{\infty}, \\
+ \Delta &= \|x^{cpu}-x^{gpu}\|_{\infty},
+ \end{aligned}
\end{equation}
where $x^{cpu}$ and $x^{gpu}$ are the solutions computed, respectively, on the CPU cluster and on the GPU cluster.
We can see that the precision of the solutions computed on the GPU cluster are sufficient, they are about $10^{-10}$,
in most numerical problems. This generator uses the sparse matrices of the Davis collection as the initial
matrices to build the large band matrices. It is executed in parallel by all the MPI processes of the cluster
so that each process constructs its own sub-matrix as a rectangular block of the global sparse matrix. Each process
-$i$ computes the size $n_i$ and the offset $offset_i$ of its sub-matrix in the global sparse matrix according to the
+$i$ computes the size $n_i$ and the offset $\Offset_i$ of its sub-matrix in the global sparse matrix according to the
size $n$ of the linear system to be solved and the number of the GPU computing nodes $p$ as follows:
\begin{equation}
n_i = \frac{n}{p},
\end{equation}
\begin{equation}
- offset_i = \left\{
- \begin{array}{l}
- 0\mbox{~if~}i=0,\\
- offset_{i-1}+n_{i-1}\mbox{~otherwise.}
- \end{array}
- \right.
+ \Offset_i =
+ \begin{cases}
+ 0 & \text{if $i=0$,}\\
+ \Offset_{i-1}+n_{i-1} & \text{otherwise.}
+ \end{cases}
\end{equation}
So each process $i$ performs several copies of the same initial matrix chosen from the Davis collection and it
puts all these copies on the main diagonal of the global matrix in order to construct a band matrix. Moreover,
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
-Matrix & $Time_{cpu}$ & $Time_{gpu}$ & $\tau$ & \# $iter$& $prec$ & $\Delta$ \\ \hline \hline
+Matrix & $\Time{cpu}$ & $\Time{gpu}$ & $\tau$ & \# $\Iter$& $\Prec$ & $\Delta$ \\ \hline \hline
2cubes\_sphere & 3.683 s & 0.870 s & 4.23 & 21 & 2.11e-14 & 8.67e-18 \\
ecology2 & 2.570 s & 0.424 s & 6.06 & 21 & 4.88e-13 & 2.08e-14 \\
finan512 & 2.727 s & 0.533 s & 5.11 & 17 & 3.22e-12 & 8.82e-14 \\
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
-Matrix & $Time_{cpu}$ & $Time_{gpu}$ & $\tau$& \# $iter$& $prec$ & $\Delta$ \\ \hline \hline
+Matrix & $\Time{cpu}$ & $\Time{gpu}$ & $\tau$& \# $\Iter$& $\Prec$ & $\Delta$ \\ \hline \hline
2cubes\_sphere & 3.597 s & 0.514 s & 6.99 & 21 & 2.11e-14 & 8.67e-18 \\
ecology2 & 2.549 s & 0.288 s & 8.83 & 21 & 4.88e-13 & 2.08e-14 \\
finan512 & 2.660 s & 0.377 s & 7.05 & 17 & 3.22e-12 & 8.82e-14 \\
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
-Matrix & $Time_{cpu}$ & $Time_{gpu}$ & $\tau$ & \# $iter$& $prec$ & $\Delta$ \\ \hline \hline
+Matrix & $\Time{cpu}$ & $\Time{gpu}$ & $\tau$ & \# $\Iter$& $\Prec$ & $\Delta$ \\ \hline \hline
2cubes\_sphere & 15.963 s & 7.250 s & 2.20 & 58 & 6.23e-16 & 3.25e-19 \\
ecology2 & 3.549 s & 2.176 s & 1.63 & 21 & 4.78e-15 & 1.06e-15 \\
finan512 & 3.862 s & 1.934 s & 1.99 & 17 & 3.21e-14 & 8.43e-17 \\
and $\lambda_j$ are, respectively, the cost and the connectivity of the cut hyperedge $e_j$. In addition,
the hypergraph partitioning is constrained to maintain the load balance between the $K$ parts:
\begin{equation}
-W_k\leq (1+\epsilon)W_{avg}\mbox{,~}(1\leq k\leq K)\mbox{~and~}(0<\epsilon<1),
+W_k\leq (1+\epsilon)\Wavg\mbox{,~}(1\leq k\leq K)\mbox{~and~}(0<\epsilon<1),
\end{equation}
-where $W_k$ is the sum of the vertex weights in the subset $\mathcal{V}_k$, $W_{avg}$ is the average part's
+where $W_k$ is the sum of the vertex weights in the subset $\mathcal{V}_k$, $\Wavg$ is the average part's
weight and $\epsilon$ is the maximum allowed imbalanced ratio.
The hypergraph partitioning is a NP-complete problem but software tools using heuristics are developed, for
Zoltan tool~\cite{ref20,ref21}. The parameters of the hypergraph partitioning are initialized as follows:
\begin{itemize}
\item The weight $w_i$ of each vertex $v_i$ is set to the number of the nonzero values on the matrix row $i$,
-\item For simplicity sake, the cost $c_j$ of each hyperedge $e_j$ is set to 1,
+\item For simplicity's sake, the cost $c_j$ of each hyperedge $e_j$ is set to 1,
\item The maximum imbalanced ratio $\epsilon$ is limited to 10\%.
\end{itemize}
We can notice from Table~\ref{tab:08} that the execution times on the cluster of 12 GPUs are significantly
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
-Matrix & $Time_{cpu}$ & $Time_{gpu}$ & $\tau$ & \# iter & $prec$ & $\Delta$ \\ \hline \hline
+Matrix & $\Time{cpu}$ & $\Time{gpu}$ & $\tau$ & \# $\Iter$ & $\Prec$ & $\Delta$ \\ \hline \hline
2cubes\_sphere & 16.430 s & 2.840 s & 5.78 & 58 & 6.23e-16 & 3.25e-19 \\
ecology2 & 3.152 s & 0.367 s & 8.59 & 21 & 4.78e-15 & 1.06e-15 \\
finan512 & 3.672 s & 0.723 s & 5.08 & 17 & 3.21e-14 & 8.43e-17 \\
\end{center}
\end{table}
-Table~\ref{tab:09} shows in the second, third and fourth columns the total communication volume on a cluster of 12 GPUs by using row-by-row partitioning or hypergraph partitioning and compressed format. The total communication volume defines the total number of the vector elements exchanged between the 12 GPUs. From these columns we can see that the two heuristics, compressed format for the vectors and the hypergraph partitioning, minimize the number of vector elements to be exchanged over the GPU cluster. The compressed format allows the GPUs to exchange the needed vector elements witout any communication overheads. The hypergraph partitioning allows to split the large sparse matrices so as to minimize data dependencies between the GPU computing nodes. However, we can notice in the fifth column that the hypergraph partitioning takes longer than the computation times. As we have mentioned before, the hypergraph partitioning method is less efficient in terms of memory consumption and partitioning time than its graph counterpart. So for the applications which often use the same sparse matrices, we can perform the hypergraph partitioning only once and, then, we save the traces in files to be reused several times. Therefore, this allows us to avoid the partitioning of the sparse matrices at each resolution of the linear systems.
+Table~\ref{tab:09} shows in the second, third and fourth columns the total communication volume on a cluster of 12 GPUs by using row-by-row partitioning or hypergraph partitioning and compressed format. The total communication volume defines the total number of the vector elements exchanged between the 12 GPUs. From these columns we can see that the two heuristics, compressed format for the vectors and the hypergraph partitioning, minimize the number of vector elements to be exchanged over the GPU cluster. The compressed format allows the GPUs to exchange the needed vector elements without any communication overheads. The hypergraph partitioning allows to split the large sparse matrices so as to minimize data dependencies between the GPU computing nodes. However, we can notice in the fifth column that the hypergraph partitioning takes longer than the computation times. As we have mentioned before, the hypergraph partitioning method is less efficient in terms of memory consumption and partitioning time than its graph counterpart. So for the applications which often use the same sparse matrices, we can perform the hypergraph partitioning only once and, then, we save the traces in files to be reused several times. Therefore, this allows us to avoid the partitioning of the sparse matrices at each resolution of the linear systems.
\begin{table}
\begin{center}
\begin{tabular}{|c||c|c|c||c|c|c|}
\hline
\multirow{2}{*}{Matrix} & \multicolumn{3}{c||}{GPU version} & \multicolumn{3}{c|}{CPU version} \\ \cline{2-7}
- & $Time_{comput}$ & $Time_{comm}$ & $Ratio$ & $Time_{comput}$ & $Time_{comm}$ & $Ratio$ \\ \hline \hline
+ & $\Time{comput}$ & $\Time{comm}$ & $\Ratio$ & $\Time{comput}$ & $\Time{comm}$ & $\Ratio$ \\ \hline \hline
2cubes\_sphere & 37.067 s & 1434.512 s & {\bf 0.026} & 312.061 s & 3453.931 s & {\bf 0.090}\\
ecology2 & 4.116 s & 501.327 s & {\bf 0.008} & 60.776 s & 1216.607 s & {\bf 0.050}\\
finan512 & 7.170 s & 386.742 s & {\bf 0.019} & 72.464 s & 932.538 s & {\bf 0.078}\\
\begin{tabular}{|c||c|c|c||c|c|c|}
\hline
\multirow{2}{*}{Matrix} & \multicolumn{3}{c||}{GPU version} & \multicolumn{3}{c|}{CPU version} \\ \cline{2-7}
- & $Time_{comput}$ & $Time_{comm}$ & $Ratio$ & $Time_{comput}$ & $Time_{comm}$ & $Ratio$ \\ \hline \hline
+ & $\Time{comput}$ & $\Time{comm}$ & $\Ratio$ & $\Time{comput}$ & $\Time{comm}$ & $\Ratio$ \\ \hline \hline
2cubes\_sphere & 27.386 s & 154.861 s & {\bf 0.177} & 342.255 s & 42.100 s & {\bf 8.130}\\
ecology2 & 3.822 s & 53.131 s & {\bf 0.072} & 69.956 s & 15.019 s & {\bf 4.658}\\
finan512 & 6.366 s & 41.155 s & {\bf 0.155} & 79.592 s & 8.604 s & {\bf 9.251}\\
\begin{tabular}{|c||c|c|c||c|c|c|}
\hline
\multirow{2}{*}{Matrix} & \multicolumn{3}{c||}{GPU version} & \multicolumn{3}{c|}{CPU version} \\ \cline{2-7}
- & $Time_{comput}$ & $Time_{comm}$ & $Ratio$ & $Time_{comput}$ & $Time_{comm}$ & $Ratio$ \\ \hline \hline
+ & $\Time{comput}$ & $\Time{comm}$ & $\Ratio$ & $\Time{comput}$ & $\Time{comm}$ & $\Ratio$ \\ \hline \hline
2cubes\_sphere & 28.440 s & 7.768 s & {\bf 3.661} & 327.109 s & 63.788 s & {\bf 5.128}\\
ecology2 & 3.652 s & 0.757 s & {\bf 4.823} & 63.632 s & 13.520 s & {\bf 4.707}\\
finan512 & 7.579 s & 4.569 s & {\bf 1.659} & 74.120 s & 22.505 s & {\bf 3.294}\\
Figure~\ref{fig:09} presents the weak scaling of four versions of the parallel GMRES algorithm on a GPU cluster. We fixed the size of a sub-matrix to 5 million of rows per GPU computing node. We used matrices having five bands generated from the symmetric matrix thermal2. This figure shows that the parallel GMRES algorithm, in its naive version or using either the compression format for vectors or the hypergraph partitioning, is not scalable on a GPU cluster due to the large amount of communications between GPUs. In contrast, we can see that the algorithm using both optimization techniques is fairly scalable. That means that in this version the cost of communications is relatively constant regardless the number of computing nodes in the cluster.\\
- Finally, as far as we are concerned, the parallel solving of a linear system can be easy to optimize when the associated matrix is regular. This is unfortunately not the case for many real-world applications. When the matrix has an irregular structure, the amount of communication between processors is not the same. Another important parameter is the size of the matrix bandwidth which has a huge influence on the amount of communication. In this work, we have generated different kinds of matrices in order to analyze several difficulties. With a bandwidth as large as possible, involving communications between all processors, which is the most difficult situation, we proposed to use two heuristics. Unfortunately, there is no fast method that optimizes the communication in any situation. For systems of non linear equations, there are different algorithms but most of them consist in linearizing the system of equations. In this case, a linear system needs to be solved. The big interest is that the matrix is the same at each step of the non linear system solving, so the partitioning method which is a time consuming step is performed only once.
+ Finally, as far as we are concerned, the parallel solving of a linear system can be easy to optimize when the associated matrix is regular. This is unfortunately not the case for many real-world applications. When the matrix has an irregular structure, the amount of communications between processors is not the same. Another important parameter is the size of the matrix bandwidth which has a huge influence on the amount of communications. In this work, we have generated different kinds of matrices in order to analyze several difficulties. With a bandwidth as large as possible, involving communications between all processors, which is the most difficult situation, we proposed to use two heuristics. Unfortunately, there is no fast method that optimizes the communications in any situation. For systems of non linear equations, there are different algorithms but most of them consist in linearizing the system of equations. In this case, a linear system needs to be solved. The big interest is that the matrix is the same at each step of the non linear system solving, so the partitioning method which is a time consuming step is performed only once.
-Another very important issue, which might be ignored by too many people, is that the communications have a greater influence on a cluster of GPUs than on a cluster of CPUs. There are two reasons for that. The first one comes from the fact that with a cluster of GPUs, the CPU/GPU data transfers slow down communications between two GPUs that are not on the same machines. The second one is due to the fact that with GPUs the ratio of the computation time over the communication time decreases since the computation time is reduced. So the impact of the communications between GPUs might be a very important issue that can limit the scalability of a parallel algorithm.
+Another very important issue, which might be ignored by too many people, is that the communications have a greater influence on a cluster of GPUs than on a cluster of CPUs. There are two reasons for that. The first one comes from the fact that with a cluster of GPUs, the CPU/GPU data transfers slow down communications between two GPUs that are not on the same machines. The second one is due to the fact that with GPUs the ratio of the computation time over the communication time decreases since the computation time is reduced. So the impact of the communications between GPUs might be a very important issue that can limit the scalability of a parallel algorithms.
%%--------------------%%
%% SECTION 7 %%
