\usepackage{url}
\usepackage{mdwlist}
\usepackage{multirow}
+\usepackage{color}
\date{}
In this paper, we use the GMRES iterative method for solving large sparse linear systems on a cluster
of GPUs. The parallel algorithm of this method is implemented using the CUDA programming language for
the GPUs and the MPI parallel environment to distribute the computations between the different GPU nodes
-of the cluster. Particularly, we focus on improving the performances of the parallel sparse matrix-vector
-multiplication. Indeed, this operation is not only very time-consuming but it also requires communications
+of the cluster. Particularly, we focus on improving the performances of the parallel sparse matrix-vector multiplication.
+Indeed, this operation is not only very time-consuming but it also requires communications
between the GPU nodes. These communications are needed to build the global vector involved in
the parallel sparse matrix-vector multiplication. It should be noted that a communication between two
GPU nodes involves data transfers between the GPU and CPU memories in the same node and the MPI communications
volume. In addition, the performances of the parallel FEM algorithm are improved by overlapping
the communication with computation.
+%%% MODIF %%%
+\textcolor{red}{ \bf Our main contribution in this work is to show the difficulties to implement the GMRES method for solving sparse linear systems on a cluster of GPUs. First, we show the main key points of the parallel GMRES algorithm on a GPU cluster. Then, we discuss the improvements of the algorithm which are mainly performed on the sparse matrix-vector multiplication when the matrix is distributed on several GPUs. In fact, on a cluster of GPUs the influence of the communications is greater than on clusters of CPUs due to the CPU/GPU communications between two GPUs that are not on the same machines. We propose to perform a hypergraph partitioning on the problem to be solved, then we reorder the matrix columns according to the partitioning scheme, and we use a compressed format for storing the vectors in such a way to minimize the communication overheads between two GPUs.}
+%%% END %%%
%%--------------------%%
%% SECTION 3 %%
%%--------------------%%
\section{{GMRES} method}
\label{sec:04}
-The generalized minimal residual method (GMRES) is an iterative method designed by Saad
-and Schultz in 1986~\cite{Saa86}. It is a generalization of the minimal residual method
-(MNRES)~\cite{Pai75} to deal with asymmetric and non Hermitian problems and indefinite
-symmetric problems.
+%%% MODIF %%%
+The generalized minimal residual method (GMRES) is an iterative method designed by Saad and Schultz in 1986~\cite{Saa86}. It is a generalization of the minimal residual method (MNRES)~\cite{Pai75} to deal with asymmetric and non Hermitian problems and indefinite symmetric problems.
Let us consider the following sparse linear system of $n$ equations:
\begin{equation}
Ax = b,
\label{eq:01}
\end{equation}
-where $A\in\mathbb{R}^{n\times n}$ is a sparse square and nonsingular matrix, $x\in\mathbb{R}^n$
-is the solution vector and $b\in\mathbb{R}^n$ is the right-hand side vector. The main idea of
-the GMRES method is to find a sequence of solutions $\{x_k\}_{k\in\mathbb{N}}$ which minimizes at
-best the residual $r_k=b-Ax_k$. The solution $x_k$ is computed in a Krylov sub-space $\mathcal{K}_k(A,v_1)$:
+where $A\in\mathbb{R}^{n\times n}$ is a sparse square and nonsingular matrix, $x\in\mathbb{R}^n$ is the solution vector and $b\in\mathbb{R}^n$ is the right-hand side vector. The main idea of the GMRES method is to find a sequence of solutions $\{x_k\}_{k\in\mathbb{N}}$ which minimizes at best the residual $r_k=b-Ax_k$. The solution $x_k$ is computed in a Krylov sub-space $\mathcal{K}_k(A,v_1)$:
\begin{equation}
\begin{array}{ll}
\mathcal{K}_{k}(A,v_{1}) \equiv \text{span}\{v_{1}, Av_{1}, A^{2}v_{1},..., A^{k-1}v_{1}\}, & v_{1}=\frac{r_{0}}{\|r_{0}\|_{2}},
\begin{equation}
r_{k} \perp A\mathcal{K}_{k}(A, v_{1}).
\end{equation}
-The GMRES method uses the Arnoldi method~\cite{Arn51} to build an orthonormal basis $V_k$
-and a Hessenberg upper matrix $\bar{H}_k$ of order $(k+1)\times k$:
-\begin{equation}
-\begin{array}{lll}
- V_{k}=\{v_{1},v_{2},...,v_{k}\}, & \text{such that} & \forall k>1, v_{k}=A^{k-1}v_{1},
-\end{array}
-\end{equation}
-and
-\begin{equation}
- AV_{k} = V_{k+1}\bar{H}_{k}.
-\label{eq:02}
-\end{equation}
-Then at each iteration $k$, a solution $x_k$ is computed in the Krylov sub-space $\mathcal{K}_{k}$
-spanned by $V_k$ as follows:
-\begin{equation}
-\begin{array}{ll}
- x_{k} = x_{0} + V_{k}y, & y\in\mathbb{R}^{k}.
-\end{array}
-\label{eq:03}
-\end{equation}
-From both equations~(\ref{eq:02}) and~(\ref{eq:03}), we can deduce:
-\begin{equation}
-\begin{array}{lll}
- r_{k} & = & b - A (x_{0} + V_{k}y) \\
- & = & r_{0} - AV_{k}y \\
- & = & \beta v_{1} - V_{k+1}\bar{H}_{k}y \\
- & = & V_{k+1}(\beta e_{1} - \bar{H}_{k}y),
-\end{array}
-\end{equation}
-where $\beta=\|r_{0}\|_{2}$ and $e_{1}=(1,0,...,0)$ is the first vector of the canonical basis of
-$\mathbb{R}^{k}$. The vector $y$ is computed in $\mathbb{R}^{k}$ so as to minimize at best the Euclidean
-norm of the residual $r_{k}$. This means that the following least-squares problem of size $k$ must
-be solved:
-\begin{equation}
- \underset{y\in\mathbb{R}^{k}}{min}\|r_{k}\|_{2}=\underset{y\in\mathbb{R}^{k}}{min}\|\beta e_{1}-\bar{H}_{k}y\|_{2}.
-\end{equation}
-The solution of this problem is computed thanks to the QR factorization of the Hessenberg matrix
-$\bar{H}_{k}$ by using the Givens rotations~\cite{Saa03,Saa86} such that:
-\begin{equation}
-\begin{array}{lll}
- \bar{H}_{k}=Q_{k}R_{k}, & Q_{k}\in\mathbb{R}^{(k+1)\times (k+1)}, & R_{k}\in\mathbb{R}^{(k+1)\times k},
-\end{array}
-\end{equation}
-where $Q_{k}Q_{k}^{T}=I_{k}$ and $R_{k}$ is a upper triangular matrix.
-
-The GMRES method finds a solution at most after $n$ iterations. So, the Krylov sub-space dimension
-can increase up-to the large size $n$ of the linear system. Then in order to avoid a large storage
-of the orthonormal basis $V_{k}$, the Arnoldi process is restricted at $m\ll n$ iterations and restarted
-with the last iterate $x_{m}$ as an initial guess to the next iteration. Moreover, GMRES sometimes
-does not converge or takes too many iterations to satisfy the convergence criteria. Therefore, in
-most cases, GMRES must contain a preconditioning step to improve its convergence. The preconditioning
-technique replaces the system~(\ref{eq:01}) with the modified systems:
-\begin{equation}M^{-1}Ax=M^{-1}b.\end{equation}
-
-Algorithm~\ref{alg:01} illustrates the main key points of the GMRES method with restarts.
-The linear system to be solved in this algorithm is left-preconditioned where $M$ is the preconditioning
-matrix. At each iteration $k$, the Arnoldi process is used (from line~$7$ to line~$17$ of algorithm~\ref{alg:01})
-to construct an orthonormal basis $V_m$ and a Hessenberg matrix $\bar{H}_m$ of order $(m+1)\times m$.
-Then, the least-squares problem is solved (line~$18$) to find the vector $y\in\mathbb{R}^m$ which minimizes
-the residual. Finally, the solution $x_m$ is computed in the Krylov sub-space spanned by $V_m$ (line~$19$).
-In practice, the GMRES algorithm stops when the Euclidean norm of the residual is small enough and/or
-the maximum number of iterations is reached.
+
+Algorithm~\ref{alg:01} illustrates the main key points of the GMRES method with restarts. The linear system to be solved in this algorithm is left-preconditioned where $M$ is the preconditioning matrix. The Arnoldi process~\cite{Arn51} is used (from line~$7$ to line~$17$ of algorithm~\ref{alg:01}) to construct an orthonormal basis $V_m$ and a Hessenberg matrix $\bar{H}_m$ of order $(m+1)\times m$ such that $m\ll n$. Then, the least-squares problem is solved (line~$18$) to find the vector $y\in\mathbb{R}^m$ which minimizes the residual. Finally, the solution $x_m$ is computed in the Krylov sub-space spanned by $V_m$ (line~$19$). In practice, the GMRES algorithm stops when the Euclidean norm of the residual is small enough and/or the maximum number of iterations is reached.
+%%% END %%%
\begin{algorithm}[!h]
\SetAlgoLined
initialized to $1$. For simplicity 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.
+the size of a thread-block in GPUs to $512$ threads.
+%%% MODIF %%%
+\textcolor{red}{\bf It would be noted that the same optimizations done on the GPU version of the parallel GMRES algorithm are performed on the CPU version.}
+%%% END %%%
\begin{table}[!h]
\begin{center}
\end{center}
\end{table}
-Table~\ref{tab:09} shows in the second and third columns the total communication volume on a cluster
-of 12 GPUs by using row-by-row partitioning and hypergraph partitioning, respectively. The total communication
-volume defines the total number of the vector elements exchanged between the 12 GPUs. This table shows
-that the hypergraph partitioning allows to split the large sparse matrices so as to minimize data
+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 the vector elements to be exchanged over the GPU cluster. The compressed format applied
+
+
+
+This table shows that 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 fourth 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.
\begin{table}[!h]
\begin{center}
-\begin{tabular}{|c|c|c|c|}
+\begin{tabular}{|c|c|c|c|c|}
\hline
-\multirow{4}{*}{Matrix} & Total communication & Total communication & Time of \\
- & volume using & volume using & hypergraph \\
- & row-by-row & hypergraph & partitioning \\
- & partitioning & partitioning & in minutes \\ \hline \hline
-2cubes\_sphere & 25 360 543 & 240 679 & 68.98 \\
-ecology2 & 26 044 002 & 73 021 & 4.92 \\
-finan512 & 26 087 431 & 900 729 & 33.72 \\
-G3\_circuit & 31 912 003 & 5 366 774 & 11.63 \\
-shallow\_water2 & 25 105 108 & 60 899 & 5.06 \\
-thermal2 & 30 012 846 & 1 077 921 & 17.88 \\ \hline \hline
-cage13 & 28 254 282 & 3 845 440 & 196.45 \\
-crashbasis & 29 020 060 & 2 401 876 & 33.39 \\
-FEM\_3D\_thermal2 & 25 263 767 & 250 105 & 49.89 \\
-language & 27 291 486 & 1 537 835 & 9.07 \\
-poli\_large & 25 053 554 & 7 388 883 & 5.92 \\
-torso3 & 25 682 514 & 613 250 & 61.51 \\ \hline
+\multirow{3}{*}{Matrix} & Total comm. vol. & Total comm. vol. & Total comm. vol. & Time of hypergraph \\
+ & using row-by row & using compressed & using hypergraph partitioning & partitioning \\
+ & partitioning & format & and compressed format & in minutes \\ \hline \hline
+2cubes\_sphere & 182 061 791 & 25 360 543 & 240 679 & 68.98 \\
+ecology2 & 181 267 000 & 26 044 002 & 73 021 & 4.92 \\
+finan512 & 182 090 692 & 26 087 431 & 900 729 & 33.72 \\
+G3\_circuit & 192 244 835 & 31 912 003 & 5 366 774 & 11.63 \\
+shallow\_water2 & 181 729 606 & 25 105 108 & 60 899 & 5.06 \\
+thermal2 & 191 350 306 & 30 012 846 & 1 077 921 & 17.88 \\ \hline \hline
+cage13 & 183 970 606 & 28 254 282 & 3 845 440 & 196.45 \\
+crashbasis & 182 931 818 & 29 020 060 & 2 401 876 & 33.39 \\
+FEM\_3D\_thermal2 & 182 503 894 & 25 263 767 & 250 105 & 49.89 \\
+language & 183 055 017 & 27 291 486 & 1 537 835 & 9.07 \\
+poli\_large & 181 381 470 & 25 053 554 & 7 388 883 & 5.92 \\
+torso3 & 183 863 292 & 25 682 514 & 613 250 & 61.51 \\ \hline
\end{tabular}
-\caption{Total communication volume on a cluster of 12 GPUs using row-by-row or hypergraph partitioning methods}
+\caption{Total communication volume on a cluster of 12 GPUs using row-by-row or hypergraph partitioning methods and compressed vectors. The total communication volume is defined as the total number of vector elements exchanged between all GPUs of the cluster.}
\label{tab:09}
\end{center}
\end{table}
+%%% MODIF %%%
+\textcolor{red}{\bf Finally, the parallel solving of a linear system can be easy to optimize when the associated matrix is regular. This is unfortunately not the case of 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 communications. In this work, we have generated different kinds of matrices in order to analyze different difficulties. With as a large bandwidth 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 once only.
+}
+
+\textcolor{red}{
+Another very important issue is that the communications have a greater influence on a cluster of GPUs than on a cluster of CPUs. There are two reasons for this. 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.}
+%%% END %%%
%%--------------------%%
%% SECTION 7 %%
