new

[GMRES_For_Journal.git] / GMRES_Journal.tex

\documentclass[11pt]{article}
%\documentclass{acmconf}

\usepackage[paper=a4paper,dvips,top=1.5cm,left=1.5cm,right=1.5cm,foot=1cm,bottom=1.5cm]{geometry}
\usepackage{times}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{amsopn}
\usepackage{xspace}
\usepackage{array}
\usepackage{epsfig}
\usepackage{algorithmic}
\usepackage[ruled,english,boxed,linesnumbered]{algorithm2e}
\usepackage[english]{algorithm2e}
\usepackage{url}
\usepackage{mdwlist}
\usepackage{multirow}

\date{}

\title{Parallel sparse linear solver with GMRES method using minimization techniques of communications for GPU clusters}

\author{
\textsc{Jacques M. Bahi}
\qquad
\textsc{Rapha\"el Couturier}\thanks{Contact author}
\qquad
\textsc{Lilia Ziane Khodja}
\mbox{}\\ %
\\
FEMTO-ST Institute, University of Franche-Comte\\
IUT Belfort-Montb\'eliard\\
Rue Engel Gros, BP 527, 90016 Belfort, \underline{France}\\
\mbox{}\\ %
\normalsize
\{\texttt{jacques.bahi},~\texttt{raphael.couturier},~\texttt{lilia.ziane\_khoja}\}\texttt{@univ-fcomte.fr}
}

\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
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
cluster.  
\end{abstract}



%%--------------------%%
%%      SECTION 1     %%
%%--------------------%%
\section{Introduction}
\label{sec:01}
Large sparse linear systems arise in  most numerical scientific or industrial simulations.
They model numerous complex problems in different areas of applications such as mathematics,
engineering, biology or physics~\cite{ref18}. However, solving these systems of equations is
often an expensive operation in terms of execution time and memory space consumption. Indeed,
the linear systems arising in most applications are very large  and have many zero
coefficients, and this sparse nature leads to irregular accesses to load the nonzero coefficients
from the memory.  

Parallel computing has become a key issue for solving sparse linear systems of large sizes.
This is due to the computing power and the storage capacity of the current parallel computers as
well as the availability of different parallel programming languages ​​and environments such as the
MPI communication standard. Nowadays, graphics processing units (GPUs) are the most commonly used
hardware accelerators in high performance computing. They are equipped with a massively parallel
architecture allowing them to compute faster than CPUs. However, the parallel computers equipped 
with GPUs introduce new programming difficulties to adapt  parallel algorithms to their architectures.

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
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
between the CPUs of the GPU nodes. For performance purposes, we propose to use a compressed storage
format to reduce the size of the vectors to be exchanged between the GPU nodes and a hypergraph partitioning
of the sparse matrix to reduce the total communication volume. 

The present paper is organized as follows. In Section~\ref{sec:02} some previous works about solving 
sparse linear systems on GPUs are presented. In Section~\ref{sec:03} is given a general overview of the GPU architectures,
followed by that the GMRES method in Section~\ref{sec:04}. In Section~\ref{sec:05} the main key points
of the parallel implementation of the GMRES method on a GPU cluster are described. Finally, in Section~\ref{sec:06}
is presented the performance improvements of the parallel GMRES algorithm on a GPU cluster.


%%--------------------%%
%%      SECTION 2     %%
%%--------------------%%
\section{Related work}
\label{sec:02}
Numerous works have shown the efficiency of GPUs for solving sparse linear systems compared
to their CPUs counterpart. Different iterative methods are implemented on one GPU, for example
Jacobi and Gauss-Seidel in~\cite{refa}, conjugate and biconjugate gradients in~\cite{refd,refe,reff,refj}
and GMRES in~\cite{refb,refc,refg,refm}. In addition, some iterative methods are implemented on 
shared memory multi-GPUs machines as~\cite{refh,refi,refk,refl}. A limited set of studies are
devoted to the parallel implementation of the iterative methods on distributed memory GPU clusters
as~\cite{refn,refo,refp}.

Traditionally, the parallel iterative algorithms do not often scale well on GPU clusters due to
the significant cost of the communications between the computing nodes. Some authors have already
studied how to reduce these communications. In~\cite{cev10}, the authors used a hypergraph partitioning
as a preprocessing to the parallel conjugate gradient algorithm in order to reduce the inter-GPU
communications over a GPU cluster. The sequential hypergraph partitioning method provided by the
PaToH tool~\cite{Cata99} is used because of the small sizes of the sparse symmetric linear systems
to be solved. In~\cite{refq}, a compression and decompression technique is proposed to reduce the
communication overheads. This technique is performed on the shared vectors to be exchanged between
the computing nodes. In~\cite{refr}, the authors studied the impact of asynchronism on parallel
iterative algorithms on local GPU clusters. Asynchronous communication primitives suppress some
synchronization barriers and allow overlap of communication and computation. In~\cite{refs}, a
communication reduction method is used for implementing finite element methods (FEM) on GPU clusters.
This method firstly uses the Reverse Cuthill-McKee reordering to reduce the total communication
volume. In addition, the performances of the parallel FEM algorithm are improved by overlapping
the communication with computation. 


%%--------------------%%
%%      SECTION 3     %%
%%--------------------%%
\section{{GPU} architectures}
\label{sec:03}
A GPU (Graphics processing unit) is a hardware accelerator for high performance computing.
Its hardware architecture is composed of hundreds of cores organized in several blocks called
\emph{streaming multiprocessors}. It is also equipped with a memory hierarchy. It has a set
of registers and a private read-write \emph{local memory} per core, a fast \emph{shared memory},
read-only \emph{constant} and \emph{texture} caches per multiprocessor and a read-write
\emph{global memory} shared by all its multiprocessors. The new architectures (Fermi, Kepler,
etc) have also L1 and L2 caches to improve the accesses to the global memory.

NVIDIA has released the CUDA platform (Compute Unified Device Architecture)~\cite{Nvi10}
which provides a high level GPGPU-based programming language (General-Purpose computing
on GPUs), allowing to program GPUs for general purpose computations. In CUDA programming
environment, all data-parallel and compute intensive portions of an application running
on the CPU are off-loaded onto the GPU. Indeed, an application developed in CUDA is a
program written in C language (or Fortran) with a minimal set of extensions to define
the parallel functions to be executed by the GPU, called \emph{kernels}. We define kernels,
as separate functions from those of the CPU, by assigning them a function type qualifiers
\verb+__global__+ or \verb+__device__+.

At the GPU level, the same kernel is executed by a large number of parallel CUDA threads
grouped together as a grid of thread blocks. Each multiprocessor of the GPU executes one
or more thread blocks in SIMD fashion (Single Instruction, Multiple Data) and in turn each
core of a GPU multiprocessor runs one or more threads within a block in SIMT fashion (Single
Instruction, Multiple threads). In order to maximize the occupation of the GPU cores, the
number of CUDA threads to be involved in a kernel execution is computed according to the
size of the problem to be solved. In contrast, the block size is restricted by the limited
memory resources of a core. On current GPUs, a thread block may contain up-to $1,024$ concurrent
threads. At any given clock cycle, the threads execute the same instruction of a kernel,
but each of them operates on different data. Moreover, threads within a block can cooperate
by sharing data through the fast shared memory and coordinate their execution through
synchronization points. In contrast, within a grid of thread blocks, there is no synchronization
at all between blocks. 

GPUs  only work on data filled in their global memory and the final results of their kernel
executions must be communicated to their hosts (CPUs). Hence, the data must be transferred
\emph{in} and \emph{out} of the GPU. However, the speed of memory copy between the CPU and
the GPU is slower than the memory copy speed of GPUs. Accordingly, it is necessary to limit
the transfer of data between the GPU and its host.


%%--------------------%%
%%      SECTION 4     %%
%%--------------------%%
\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. 

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)$:
\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}},
\end{array} 
\end{equation}
such that the Petrov-Galerkin condition is satisfied:
\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. 

\begin{algorithm}[!h]
  \SetAlgoLined
  \Entree{$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)}
  \BlankLine
  $r_{0} \leftarrow M^{-1}(b - Ax_{0})$\;
  $\beta \leftarrow \|r_{0}\|_{2}$\;
  $\alpha \leftarrow \|M^{-1}b\|_{2}$\;
  $convergence \leftarrow false$\;
  $k \leftarrow 1$\;
  \BlankLine
  \While{$(\neg convergence)$}{
    $v_{1} \leftarrow r_{0} / \beta$\;
    \For{$j=1$ {\bf to} $m$}{
      $w_{j} \leftarrow M^{-1}Av_{j}$\;
      \For{$i=1$ {\bf to} $j$}{
        $h_{i,j} \leftarrow (w_{j},v_{i})$\;
        $w_{j} \leftarrow w_{j} - h_{i,j} \times v_{i}$\;
      }
      $h_{j+1,j} \leftarrow \|w_{j}\|_{2}$\;
      $v_{j+1} \leftarrow w_{j} / h_{j+1,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}$\;
    \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$\;
    }{
      $x_{0} \leftarrow x_{m}$\;
      $r_{0} \leftarrow r_{m}$\;
      $k \leftarrow k + 1$\;
    }
  }
\caption{Left-preconditioned GMRES algorithm with restarts}
\label{alg:01}
\end{algorithm}



%%--------------------%%
%%      SECTION 5     %%
%%--------------------%%
\section{Parallel GMRES method on {GPU} clusters}
\label{sec:05}

\subsection{Parallel implementation on a GPU cluster}
\label{sec:05.01}
The implementation of the GMRES algorithm on a GPU cluster is performed by using
a parallel heterogeneous programming. We use the programming language CUDA for the
GPUs and the parallel environment MPI for the distribution of the computations between
the GPU computing nodes. In this work, a GPU computing node is composed of a GPU and
a CPU core managed by a MPI process.

Let us consider a cluster composed of $p$ GPU computing nodes. First, the sparse linear
system~(\ref{eq:01}) is split into $p$ sub-linear systems, each is attributed to a GPU
computing node. We partition row-by-row the sparse matrix $A$ and both vectors $x$ and
$b$ in $p$ parts (see Figure~\ref{fig:01}). The data issued from the partitioning operation
are off-loaded on the GPU global memories to be proceeded by the GPUs. Then, all the
computing nodes of the GPU cluster execute the same GMRES iterative algorithm but on
different data. Finally, the GPU computing nodes synchronize their computations by using
MPI communication routines to solve the global sparse linear system. In what follows,
the computing nodes sharing data are called the neighboring nodes.

\begin{figure}
\centering
  \includegraphics[width=80mm,keepaspectratio]{Figures/partition}
\caption{Data partitioning of the sparse matrix $A$, the solution vector $x$ and the right-hand side $b$ in $4$ partitions}
\label{fig:01}
\end{figure}

In order to exploit the computing power of the GPUs, we have to execute maximum computations
on GPUs to avoid the data transfers between the GPU and its host (CPU), and to maximize the
GPU cores utilization to hide global memory access latency. The implementation of the GMRES
algorithm is performed by executing the functions operating on vectors and matrices as kernels
on GPUs. These operations are often easy to parallelize and more efficient on parallel architectures
when they operate on large vectors. We use the fastest routines of the CUBLAS library
(CUDA Basic Linear Algebra Subroutines) to implement the dot product (\verb+cublasDdot()+),
the Euclidean norm (\verb+cublasDnrm2()+) and the AXPY operation (\verb+cublasDaxpy()+).
In addition, we have coded in CUDA a kernel for the scalar-vector product (lines~$7$ and~$15$
of Algorithm~\ref{alg:01}), a kernel for solving the least-squares problem (line~$18$) and a
kernel for solution vector updates (line~$19$).

The solution of the least-squares problem in the GMRES algorithm is based on:
\begin{itemize}
  \item a QR factorization of the Hessenberg matrix $\bar{H}$ by using plane rotations and,
  \item backward-substitution method to compute the vector $y$ minimizing the residual.
\end{itemize}
This operation is not easy to parallelize and it is not interesting to implement it on GPUs. 
However, the size $m$ of the linear least-squares problem to solve in the GMRES method with
restarts is very small. So, this problem is solved in sequential by one GPU thread.

The most important operation in the GMRES method is the sparse matrix-vector multiplication.
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
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
format and remaining entries of exceptional rows in COO format. It combines the efficiency of
ELL, due to the regularity of its memory accessing and the flexibility of COO which is insensitive
to the matrix structure. 

In the parallel GMRES algorithm, the GPU computing nodes must exchange between them their shared data in
order to construct the global vector necessary to compute the parallel sparse matrix-vector
multiplication (SpMV). In fact, each computing node has locally the vector elements corresponding
to the rows of its sparse sub-matrix and, in order to compute its part of the SpMV, it also
requires the vector elements of its neighboring nodes corresponding to the column indices in
which its local sub-matrix has nonzero values. Consequently, each computing node manages a global
vector composed of a local vector of size $\frac{n}{p}$ and a shared vector of size $S$:
\begin{equation}
  S = bw - \frac{n}{p},
\label{eq:11}
\end{equation}
where $\frac{n}{p}$ is the size of the local vector and $bw$ is the bandwidth of the local sparse
sub-matrix which represents the number of columns between the minimum and the maximum column indices
(see Figure~\ref{fig:01}). In order to improve memory accesses, we use the texture memory to
cache elements of the global vector.

On a GPU cluster, the exchanges of the shared vectors elements between the neighboring nodes are
performed as follows:
\begin{itemize}
  \item at the level of the sending node: data transfers of the shared data from the GPU global memory
to the CPU memory by using the CUBLAS communication routine \verb+cublasGetVector()+,
  \item data exchanges between the CPUs by the MPI communication routine \verb+MPI_Alltoallv()+ and,
  \item at the level of the receiving node: data transfers of the received shared data from the CPU
memory to the GPU global memory by using CUBLAS communication routine \verb+cublasSetVector()+. 
\end{itemize}  

\subsection{Experimentations}
\label{sec:05.02}
The experiments are done on a cluster composed of six machines interconnected by an Infiniband network
of $20$~GB/s. Each machine is a Xeon E5530 Quad-Core running at $2.4$~GHz. It provides $12$~GB of RAM
memory with a memory bandwidth of $25.6$~GB/s and it is equipped with two Tesla C1060 GPUs. Each GPU
is composed of $240$ cores running at $1.3$ GHz and has $4$~GB of global memory with a memory bandwidth
of $102$~GB/s. The GPU is connected to the CPU via a PCI-Express 16x Gen2.0 with a throughput of $8$~GB/s.
Figure~\ref{fig:02} shows the general scheme of the GPU cluster.

\begin{figure}[!h]
\centering
  \includegraphics[width=80mm,keepaspectratio]{Figures/clusterGPU}
\caption{A cluster composed of six machines, each equipped with two Tesla C1060 GPUs}
\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
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
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$).

\begin{table}
\begin{center}
\begin{tabular}{|c|c|r|r|r|} 
\hline
Matrix type                 & Name              & \# Rows   & \# Nonzeros  & Bandwidth \\\hline \hline
\multirow{6}{*}{Symmetric}  & 2cubes\_sphere    & 101 492   & 1 647 264    & 100 464        \\
                            & ecology2          & 999 999   & 4 995 991    & 2 001          \\ 
                            & finan512          & 74 752    & 596 992      & 74 725         \\ 
                            & G3\_circuit       & 1 585 478 & 7 660 826    & 1 219 059      \\
                            & shallow\_water2   & 81 920    & 327 680      & 58 710         \\
                            & thermal2          & 1 228 045 & 8 580 313    & 1 226 629      \\ \hline \hline
\multirow{6}{*}{Asymmetric} & cage13            & 445 315   & 7 479 343    & 318 788        \\
                            & crashbasis        & 160 000   & 1 750 416    & 120 202        \\
                            & FEM\_3D\_thermal2 & 147 900   & 3 489 300    & 117 827        \\
                            & language          & 399 130   & 1 216 334    & 398 622        \\
                            & poli\_large       & 15 575    & 33 074       & 15 575         \\
                            & torso3            & 259 156   & 4 429 042    & 216 854        \\ \hline
\end{tabular}
\caption{Main characteristics of the sparse matrices chosen from the Davis collection}
\label{tab:01}
\end{center}
\end{table}

\begin{figure}
\centering
  \includegraphics[width=120mm,keepaspectratio]{Figures/matrices}
\caption{Structures of the sparse matrices chosen from the Davis collection}
\label{fig:03}
\end{figure}

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 
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
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{table}[!h]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|} 
\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 \\
G3\_circuit       & 1.016 s     & 0.649 s      & 1.56  & 22     & 1.04e-12      & 2.00e-15 \\
shallow\_water2   & 0.061 s     & 0.044 s      & 1.38  & 17     & 5.42e-22      & 2.71e-25 \\
thermal2          & 1.666 s     & 0.880 s      & 1.89  & 21     & 6.58e-12      & 2.77e-16 \\ \hline \hline
cage13            & 0.721 s     & 0.338 s      & 2.13  & 26     & 3.37e-11      & 2.66e-15 \\
crashbasis        & 1.349 s     & 0.830 s      & 1.62  & 121    & 9.10e-12      & 6.90e-12 \\
FEM\_3D\_thermal2 & 0.797 s     & 0.419 s      & 1.90  & 64     & 3.87e-09      & 9.09e-13 \\
language          & 2.252 s     & 1.204 s      & 1.87  & 90     & 1.18e-10      & 8.00e-11 \\
poli\_large       & 0.097 s     & 0.095 s      & 1.02  & 69     & 4.98e-11      & 1.14e-12 \\
torso3            & 4.242 s     & 2.030 s      & 2.09  & 175    & 2.69e-10      & 1.78e-14 \\ \hline
\end{tabular}
\caption{Performances of the parallel GMRES algorithm on a cluster of 24 CPU cores vs. a cluster of 12 GPUs}
\label{tab:02}
\end{center}
\end{table}

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
is computed as follows:
\begin{equation}
  \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
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}
\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}$,
and the parallel GMRES algorithm computes almost the same solutions in both CPU and GPU clusters, with $\Delta$ varying
from $10^{-11}$ to $10^{-25}$.

Afterwards, we evaluate the performances of the parallel GMRES algorithm for solving large linear systems. We have
developed in C programming language a generator of large sparse matrices having a band structure which arises 
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
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.
\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,
it fulfills the empty spaces between two successive copies by small copies, \textit{lower\_copy} and \textit{upper\_copy},
of the same initial matrix. Figure~\ref{fig:04} shows a generation of a sparse bended matrix by four computing nodes.

\begin{figure}
\centering
  \includegraphics[width=100mm,keepaspectratio]{Figures/generation}
\caption{Example of the generation of a large sparse and band matrix by four computing nodes.}
\label{fig:04}
\end{figure}

Table~\ref{tab:03} shows the main characteristics (the number of nonzero values and the bandwidth) of the
large sparse matrices generated from those of the Davis collection. These matrices are associated to the
linear systems of 25 million of unknown values (each generated sparse matrix has 25 million rows). In Table~\ref{tab:04}
we show the performances of the parallel GMRES algorithm for solving large linear systems associated to the
sparse band matrices of Table~\ref{tab:03}. The fourth column gives the ratio between the execution time
spent on a cluster of 24 CPU cores and that spent on a cluster of 12 GPUs. We can notice from these ratios
that for solving large sparse matrices the GPU cluster is more efficient (about 5 times faster) than the CPU
cluster. The computing power of the GPUs allows to accelerate the computation of the functions operating
on large vectors of the parallel GMRES algorithm.

\begin{table}[!h]
\begin{center}
\begin{tabular}{|c|c|r|r|} 
\hline
Matrix type                 & Name              & \# nonzeros & Bandwidth \\ \hline \hline
\multirow{6}{*}{Symmetric}  & 2cubes\_sphere    & 413 703 602 & 198 836   \\
                            & ecology2          & 124 948 019 & 2 002     \\ 
                            & finan512          & 278 175 945 & 123 900   \\
                            & G3\_circuit       & 125 262 292 & 1 891 887 \\            
                            & shallow\_water2   & 100 235 292 & 62 806    \\
                            & thermal2          & 175 300 284 & 2 421 285 \\ \hline \hline
\multirow{6}{*}{Asymmetric} & cage13            & 435 770 480 & 352 566   \\
                            & crashbasis        & 409 291 236 & 200 203   \\
                            & FEM\_3D\_thermal2 & 595 266 787 & 206 029   \\
                            & language          & 76 912 824  & 398 626   \\
                            & poli\_large       & 53 322 580  & 15 576    \\
                            & torso3            & 433 795 264 & 328 757   \\ \hline
\end{tabular}
\caption{Main characteristics of the sparse and band matrices generated from the sparse matrices of the Davis collection.}
\label{tab:03}
\end{center}
\end{table}


\begin{table}[!h]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|} 
\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 \\
G3\_circuit       & 4.656 s     & 1.024 s      & 4.54   & 22       & 1.04e-12 & 5.00e-15 \\
shallow\_water2   & 2.268 s     & 0.384 s      & 5.91   & 17       & 5.54e-21 & 7.92e-24 \\
thermal2          & 4.650 s     & 1.130 s      & 4.11   & 21       & 8.89e-12 & 3.33e-16 \\ \hline \hline
cage13            & 6.068 s     & 1.054 s      & 5.75   & 26       & 3.29e-11 & 1.59e-14 \\
crashbasis        & 25.906 s    & 4.569 s      & 5.67   & 135      & 6.81e-11 & 4.61e-15 \\
FEM\_3D\_thermal2 & 13.555 s    & 2.654 s      & 5.11   & 64       & 3.88e-09 & 1.82e-12 \\
language          & 13.538 s    & 2.621 s      & 5.16   & 89       & 2.11e-10 & 1.60e-10 \\
poli\_large       & 8.619 s     & 1.474 s      & 5.85   & 69       & 5.05e-11 & 6.59e-12 \\
torso3            & 35.213 s    & 6.763 s      & 5.21   & 175      & 2.69e-10 & 2.66e-14 \\ \hline
\end{tabular}
\caption{Performances of the parallel GMRES algorithm for solving large sparse linear systems associated
to band matrices on a cluster of 24 CPU cores vs. a cluster of 12 GPUs.}
\label{tab:04}
\end{center}
\end{table}


%%--------------------%%
%%      SECTION 6     %%
%%--------------------%%
\section{Minimization of communications}
\label{sec:06}
The parallel sparse matrix-vector multiplication requires  data exchanges between the GPU computing nodes
to construct the global vector. However, a GPU cluster requires communications between the GPU nodes and the
data transfers between the GPUs and their hosts CPUs. In fact, a communication between two GPU nodes implies:  
a data transfer from the GPU memory to the CPU memory at the sending node, a MPI communication between the CPUs
of two GPU nodes, and a data transfer from the CPU memory to the GPU memory at the receiving node. Moreover,
the data transfers between a GPU and a CPU are considered as the most expensive communications on a GPU cluster.
For example in our GPU cluster, the data throughput between a GPU and a CPU is of 8 GB/s which is about twice
lower than the data transfer rate between CPUs (20 GB/s) and 12 times lower than the memory bandwidth of the
GPU global memory (102 GB/s). In this section, we propose two solutions to improve the execution time of the
parallel GMRES algorithm on GPU clusters.

\subsection{Compressed storage format of the sparse vectors}
\label{sec:06.01}
In Section~\ref{sec:05.01}, the SpMV multiplication uses a global vector having a size equivalent to the matrix
bandwidth (see Formula~\ref{eq:11}). However, we can notice that a SpMV multiplication does not often need all
the vector elements of the global vector composed of the local and shared sub-vectors. For example in Figure~\ref{fig:01},
 node 1 only needs  a single vector element from  node 0 (element 1), two elements from node 2 (elements 8
and 9) and two elements from  node 3 (elements 13 and 14). Therefore to reduce the communication overheads
of the unused vector elements, the GPU computing nodes must exchange between them only the vector elements necessary
to perform their local sparse matrix-vector multiplications.

\begin{figure}
\centering
  \includegraphics[width=120mm,keepaspectratio]{Figures/compress}
\caption{Example of data exchanges between node 1 and its neighbors 0, 2 and 3.}
\label{fig:05}
\end{figure}

We propose to use a compressed storage format of the sparse global vector. In Figure~\ref{fig:05}, we show an
example of the data exchanges between node 1 and its neighbors to construct the compressed global vector.
First, the neighboring nodes 0, 2 and 3 determine the vector elements needed by node 1 and, then, they send
only these elements to it. Node 1 receives these shared elements in a compressed vector. However to compute
the sparse matrix-vector multiplication, it must first copy the received elements to the corresponding indices
in the global vector. In order to avoid this process at each iteration, we propose to reorder the columns of the
local sub-matrices so as to use the shared vectors in their compressed storage format (see Figure~\ref{fig:06}).
For performance purposes, the computation of the shared data to send to the neighboring nodes is performed
by the GPU as a kernel. In addition, we use the MPI point-to-point communication routines: a blocking send routine
\verb+MPI_Send()+ and a nonblocking receive routine \verb+MPI_Irecv()+.

\begin{figure}
\centering
  \includegraphics[width=100mm,keepaspectratio]{Figures/reorder}
\caption{Reordering of the columns of a local sparse matrix.}
\label{fig:06}
\end{figure}

Table~\ref{tab:05} shows the performances of the parallel GMRES algorithm using the compressed storage format
of the sparse global vector. The results are obtained from solving large linear systems associated to the sparse
band matrices presented in Table~\ref{tab:03}. We can see from Table~\ref{tab:05} that the execution times 
of the parallel GMRES algorithm on a cluster of 12 GPUs are improved by about 38\% compared to those presented
in Table~\ref{tab:04}. In addition, the ratios between the execution times spent on the cluster of 24 CPU cores
and those spent on the cluster of 12 GPUs have increased. Indeed, the reordering of the sparse sub-matrices and
the use of a compressed storage format for the sparse vectors minimize the communication overheads between the
GPU computing nodes.

\begin{table}[!h]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|} 
\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 \\
G3\_circuit       & 3.139 s     & 0.480 s      & 6.53  & 22       & 1.04e-12 & 5.00e-15 \\
shallow\_water2   & 2.195 s     & 0.253 s      & 8.68  & 17       & 5.54e-21 & 7.92e-24 \\
thermal2          & 3.206 s     & 0.463 s      & 6.93  & 21       & 8.89e-12 & 3.33e-16 \\ \hline \hline
cage13            & 5.560 s     & 0.663 s      & 8.39  & 26       & 3.29e-11 & 1.59e-14 \\
crashbasis        & 25.802 s    & 3.511 s      & 7.35  & 135      & 6.81e-11 & 4.61e-15 \\
FEM\_3D\_thermal2 & 13.281 s    & 1.572 s      & 8.45  & 64       & 3.88e-09 & 1.82e-12 \\
language          & 12.553 s    & 1.760 s      & 7.13  & 89       & 2.11e-10 & 1.60e-10 \\
poli\_large       & 8.515 s     & 1.053 s      & 8.09  & 69       & 5.05e-11 & 6.59e-12 \\
torso3            & 31.463 s    & 3.681 s      & 8.55  & 175      & 2.69e-10 & 2.66e-14 \\ \hline
\end{tabular}
\caption{Performances of the parallel GMRES algorithm using a compressed storage format of the sparse
vectors for solving large sparse linear systems associated to band matrices on a cluster of 24 CPU cores vs.
a cluster of 12 GPUs.}
\label{tab:05}
\end{center}
\end{table}


\subsection{Hypergraph partitioning}
\label{sec:06.02}
In this section, we use another structure of the sparse matrices. We are interested in sparse matrices
whose nonzero values are distributed along their large bandwidths. We developed in C programming
language a generator of sparse matrices having five bands (see Figure~\ref{fig:07}). The principle of
this generator is the same as the one presented in Section~\ref{sec:05.02}. However, the copies made from the
initial sparse matrix, chosen from the Davis collection, are placed on the main diagonal and on two
off-diagonals on the left and right of the main diagonal. Table~\ref{tab:06} shows the main characteristics
of sparse matrices of size 25 million of rows and generated from those of the Davis collection. We can
see in the fourth column that the bandwidths of these matrices are as large as their sizes.

\begin{figure}
\centering
  \includegraphics[width=100mm,keepaspectratio]{Figures/generation_1}
\caption{Example of the generation of a large sparse matrix having five bands by four computing nodes.}
\label{fig:07}
\end{figure}


\begin{table}[!h]
\begin{center}
\begin{tabular}{|c|c|r|r|} 
\hline
Matrix type                 & Name              & \# nonzeros   & Bandwidth \\ \hline \hline
\multirow{6}{*}{Symmetric}  & 2cubes\_sphere    & 829 082 728   & 24 999 999     \\
                            & ecology2          & 254 892 056   & 25 000 000     \\ 
                            & finan512          & 556 982 339   & 24 999 973     \\ 
                            & G3\_circuit       & 257 982 646   & 25 000 000     \\
                            & shallow\_water2   & 200 798 268   & 25 000 000     \\
                            & thermal2          & 359 340 179   & 24 999 998     \\ \hline \hline
\multirow{6}{*}{Asymmetric} & cage13            & 879 063 379   & 24 999 998     \\
                            & crashbasis        & 820 373 286   & 24 999 803     \\
                            & FEM\_3D\_thermal2 & 1 194 012 703 & 24 999 998     \\
                            & language          & 155 261 826   & 24 999 492     \\
                            & poli\_large       & 106 680 819   & 25 000 000     \\
                            & torso3            & 872 029 998   & 25 000 000     \\ \hline
\end{tabular}
\caption{Main characteristics of the sparse matrices having five band and generated from the sparse matrices of the Davis collection.}
\label{tab:06}
\end{center}
\end{table}

In Table~\ref{tab:07} we give the performances of the parallel GMRES algorithm on the CPU and GPU
clusters for solving large linear systems associated to the sparse matrices shown in Table~\ref{tab:06}.
We can notice from the ratios given in the fourth column that solving sparse linear systems associated
to matrices having large bandwidth on the GPU cluster is as fast as on the CPU cluster. This is due 
to the large total communication volume necessary to synchronize the computations over the cluster.
In fact, the naive partitioning row-by-row or column-by-column of this type of sparse matrices links
a GPU node to many neighboring nodes and produces a significant number of data dependencies between
the different GPU nodes.   

\begin{table}[!h]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|} 
\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 \\
G3\_circuit       & 4.636 s     & 2.811 s      & 1.65  & 22     & 1.08e-14 & 1.77e-16 \\
shallow\_water2   & 2.738 s     & 1.539 s      & 1.78  & 17     & 5.54e-23 & 3.82e-26 \\
thermal2          & 5.017 s     & 2.587 s      & 1.94  & 21     & 8.25e-14 & 4.34e-18 \\ \hline \hline
cage13            & 9.315 s     & 3.227 s      & 2.89  & 26     & 3.38e-13 & 2.08e-16 \\
crashbasis        & 35.980 s    & 14.770 s     & 2.43  & 127    & 1.17e-12 & 1.56e-17 \\
FEM\_3D\_thermal2 & 24.611 s    & 7.749 s      & 3.17  & 64     & 3.87e-11 & 2.84e-14 \\
language          & 16.859 s    & 9.697 s      & 1.74  & 89     & 2.17e-12 & 1.70e-12 \\
poli\_large       & 10.200 s    & 6.534 s      & 1.56  & 69     & 5.14e-13 & 1.63e-13 \\
torso3            & 49.074 s    & 19.397 s     & 2.53  & 175    & 2.69e-12 & 2.77e-16 \\ \hline
\end{tabular}
\caption{Performances of the parallel GMRES algorithm using a compressed storage format of the sparse
vectors for solving large sparse linear systems associated to matrices having five bands on a cluster
of 24 CPU cores vs. a cluster of 12 GPUs.}
\label{tab:07}
\end{center}
\end{table}

We propose to use a hypergraph partitioning method which is well adapted to numerous structures 
of sparse matrices~\cite{Cat99}. Indeed, it can well model the communications between the computing
nodes especially for the asymmetric and rectangular matrices. It gives in most cases good reductions
of the total communication volume. Nevertheless, it is more expensive in terms of execution time and
memory space consumption than the partitioning method based on graphs.  

The sparse matrix $A$ of the linear system to be solved is modelled as a hypergraph
$\mathcal{H}=(\mathcal{V},\mathcal{E})$ as follows:
\begin{itemize}
\item each matrix row $i$ ($0\leq i<n$) corresponds to a vertex $v_i\in\mathcal{V}$,
\item each matrix column $j$ ($0\leq j<n$) corresponds to a hyperedge $e_j\in\mathcal{E}$, such that:
$\forall a_{ij}$ is a nonzero value of the matrix $A$, $v_i\in pins[e_j]$,
\item $w_i$ is the weight of vertex $v_i$,
\item $c_j$ is the cost of hyperedge $e_j$.  
\end{itemize}
A $K$-way partitioning of a hypergraph $\mathcal{H}=(\mathcal{V},\mathcal{E})$ is defined as a set 
of $K$ pairwise disjoint non-empty subsets (or parts) of the vertex set $\mathcal{V}$: $\mathcal{P}=\{\mathcal{V}_1,\ldots,\mathcal{V}_k\}$,
such that $\mathcal{V}=\displaystyle\cup_{k=1}^K\mathcal{V}_{k}$. Each computing node is in charge of
a vertex subset. Figure~\ref{fig:08} shows an example of a hypergraph partitioning of a sparse matrix
of size $(9\times 9)$ into three parts. The circles and squares correspond, respectively, to the vertices
and hyperedges of the hypergraph. The solid squares define the cut hyperedges connecting at least two
different parts. The connectivity $\lambda_j$ denotes the number of different parts spanned by the cut
hyperedge $e_j$.

\begin{figure}
\centering
  \includegraphics[width=130mm,keepaspectratio]{Figures/hypergraph}
\caption{A hypergraph partitioning of a sparse matrix between three computing nodes.}
\label{fig:08}
\end{figure}

The cut hyperedges model the communications between the different GPU computing nodes in the cluster,
necessary to perform the SpMV multiplication. Indeed, each hyperedge $e_j$ defines a set of atomic
computations $b_i\leftarrow b_i+a_{ij}x_j$ of the SpMV multiplication which needs the $j^{th}$ element
of  vector $x$. Therefore pins of hyperedge $e_j$ ($pins[e_j]$) denote the set of matrix rows requiring
the same vector element $x_j$. For example in Figure~\ref{fig:08}, hyperedge $e_9$ whose pins are:
$pins[e_9]=\{v_2,v_5,v_9\}$ represents matrix rows 2, 5 and 9 requiring the vector element $x_9$
to compute in parallel the atomic operations: $b_2\leftarrow b_2+a_{29}x_9$, $b_5\leftarrow b_5+a_{59}x_9$
and $b_9\leftarrow b_9+a_{99}x_9$. However, $x_9$ is a vector element of the computing node 3 and it must
be sent to the neighboring nodes 1 and 2.

The hypergraph partitioning allows to reduce the total communication volume while maintaining the computational
load balance between the computing nodes. Indeed, it minimizes at best the following sum:
\begin{equation}
\mathcal{X}(\mathcal{P}) = \displaystyle\sum_{e_j\in\mathcal{E}_C} c_j(\lambda_j-1),
\end{equation}
where $\mathcal{E}_C$ is the set of the cut hyperedges issued from the partitioning $\mathcal{P}$, $c_j$
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),
\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
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
example: hMETIS~\cite{Kar98}, PaToH~\cite{Cata99} and Zoltan~\cite{Dev06}. Due to the large sizes of the
linear systems to be solved, we use a parallel hypergraph partitioning which must be performed by at least 
two MPI processes. The hypergraph model $\mathcal{H}$ of the sparse matrix is split into $p$ (number of computing
nodes) sub-hypergraphs $\mathcal{H}_k=(\mathcal{V}_k,\mathcal{E}_k)$, $0\leq k<p$, then the parallel partitioning
is applied by using the MPI communication routines.

Table~\ref{tab:08} shows the performances of the parallel GMRES algorithm for solving the linear systems
associated to the sparse matrices presented in Table~\ref{tab:06}. In the experiments, we have used the
compressed storage format of the sparse vectors and the parallel hypergraph partitioning developed in the
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 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
improved compared to those presented in Table~\ref{tab:07}. The hypergraph partitioning applied on the large
sparse matrices having large bandwidths have improved the execution times on the GPU cluster by about 65\%.

\begin{table}[!h]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|} 
\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 \\
G3\_circuit       & 4.468 s     & 0.971 s      & 4.60 & 22     & 1.08e-14 & 1.77e-16 \\ 
shallow\_water2   & 2.647 s     & 0.312 s      & 8.48 & 17     & 5.54e-23 & 3.82e-26 \\ 
thermal2          & 4.190 s     & 0.666 s      & 6.29 & 21     & 8.25e-14 & 4.34e-18 \\ \hline \hline
cage13            & 8.077 s     & 1.584 s      & 5.10 & 26     & 3.38e-13 & 2.08e-16 \\
crashbasis        & 35.173 s    & 5.546 s      & 6.34 & 127    & 1.17e-12 & 1.56e-17 \\
FEM\_3D\_thermal2 & 24.825 s    & 3.113 s      & 7.97 & 64     & 3.87e-11 & 2.84e-14 \\
language          & 16.706 s    & 2.522 s      & 6.62 & 89     & 2.17e-12 & 1.70e-12 \\
poli\_large       & 12.715 s    & 3.989 s      & 3.19 & 69     & 5.14e-13 & 1.63e-13 \\
torso3            & 48.459 s    & 6.234 s      & 7.77 & 175    & 2.69e-12 & 2.77e-16 \\ \hline
\end{tabular}
\caption{Performances of the parallel GMRES algorithm using a compressed storage format of the sparse
vectors and a hypergraph partitioning method for solving large sparse linear systems associated to matrices
having five bands on a cluster of 24 CPU cores vs. a cluster of 12 GPUs.}
\label{tab:08}
\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
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.
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}[!h]
\begin{center}
\begin{tabular}{|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       
\end{tabular}
\caption{Total communication volume on a cluster of 12 GPUs using row-by-row or hypergraph partitioning methods}
\label{tab:09}
\end{center}
\end{table}


%%--------------------%%
%%      SECTION 7     %%
%%--------------------%%
\section{Conclusion and perspectives}
\label{sec:07}
In this paper, we have aimed at harnessing the computing power of a GPU cluster for 
solving large sparse linear systems. We have implemented the parallel algorithm of the 
GMRES iterative method. We have used a heterogeneous parallel programming based on the
CUDA language to program the GPUs and the MPI parallel environment to distribute the
computations between the GPU nodes on the cluster.

The experiments have shown that solving large sparse linear systems is more efficient 
on a cluster of GPUs than on a cluster of CPUs. However, the efficiency of a GPU cluster
is ensured as long as the spatial and temporal localization of the data is well managed. 
The data dependency scheme on a GPU cluster is related to the sparse structures of the
matrices (positions of the nonzero values) and the number of the computing nodes. We have
shown that a large number of communications between the GPU computing nodes affects 
considerably the performances of the parallel GMRES algorithm on the GPU cluster. Therefore,
we have proposed to reorder the columns of the sparse local sub-matrices on each GPU node
and to use a compressed storage format for the sparse vector involved in the parallel
sparse matrix-vector multiplication. This solution allows to minimize the communication
overheads. In addition, we have shown that it is interesting to choose a partitioning method
according to the structure of the sparse matrix. This reduces the total communication
volume between the GPU computing nodes.

In future works, it would be interesting to implement and study the scalability of the
parallel GMRES algorithm on large GPU clusters (hundreds or thousands of GPUs) or on geographically
distant GPU clusters. In this context, other methods might be used to reduce communication
and to improve the performances of the parallel GMRES algorithm as the multisplitting methods.
The recent GPU hardware and software architectures provide the GPU-Direct system which allows
two GPUs, placed in the same machine or in two remote machines, to exchange data without using
CPUs. This improves the data transfers between GPUs. Finally, it would be interesting to implement
other iterative methods on GPU clusters for solving large sparse linear or nonlinear systems.

\paragraph{Acknowledgments}
This paper is based upon work supported by the R\'egion de Franche-Comt\'e.

\bibliography{bib}
\bibliographystyle{abbrv}
\end{document}