From c8b308b236018000b34d6d58d5bdc4cb8313111b Mon Sep 17 00:00:00 2001 From: couturie Date: Tue, 19 Mar 2013 20:59:23 +0100 Subject: [PATCH 1/1] modif lilia ch12 --- BookGPU/Chapters/chapter12/ch12.aux | 64 +- BookGPU/Chapters/chapter12/ch12.tex | 470 ++----------- BookGPU/Chapters/chapter12/ch12.tex~ | 976 ++++++++++++++++----------- BookGPU/Chapters/chapter13/ch13.aux | 154 ++--- BookGPU/Chapters/chapter16/ch16.aux | 126 ++-- BookGPU/Chapters/chapter18/ch18.aux | 70 +- 6 files changed, 802 insertions(+), 1058 deletions(-) diff --git a/BookGPU/Chapters/chapter12/ch12.aux b/BookGPU/Chapters/chapter12/ch12.aux index d3b6ca3..f25a484 100644 --- a/BookGPU/Chapters/chapter12/ch12.aux +++ b/BookGPU/Chapters/chapter12/ch12.aux @@ -52,74 +52,48 @@ \newlabel{ch12:fig:02}{{11.2}{260}} \@writefile{lof}{\contentsline {figure}{\numberline {11.3}{\ignorespaces Columns reordering of a sparse sub-matrix.\relax }}{261}} \newlabel{ch12:fig:03}{{11.3}{261}} -\@writefile{lof}{\contentsline {figure}{\numberline {11.4}{\ignorespaces General scheme of the GPU cluster of tests composed of six machines, each with two GPUs.\relax }}{262}} -\newlabel{ch12:fig:04}{{11.4}{262}} \@writefile{toc}{\contentsline {section}{\numberline {11.4}Experimental results}{262}} \newlabel{ch12:sec:04}{{11.4}{262}} -\@writefile{lof}{\contentsline {figure}{\numberline {11.5}{\ignorespaces Sketches of sparse matrices chosen from the Davis's collection.\relax }}{263}} +\@writefile{lof}{\contentsline {figure}{\numberline {11.4}{\ignorespaces General scheme of the GPU cluster of tests composed of six machines, each with two GPUs.\relax }}{262}} +\newlabel{ch12:fig:04}{{11.4}{262}} +\@writefile{lof}{\contentsline {figure}{\numberline {11.5}{\ignorespaces Sketches of sparse matrices chosen from the Davis collection.\relax }}{263}} \newlabel{ch12:fig:05}{{11.5}{263}} -\@writefile{lot}{\contentsline {table}{\numberline {11.1}{\ignorespaces Main characteristics of sparse matrices chosen from the Davis's collection.\relax }}{264}} -\newlabel{ch12:tab:01}{{11.1}{264}} +\@writefile{lot}{\contentsline {table}{\numberline {11.1}{\ignorespaces Main characteristics of sparse matrices chosen from the Davis collection.\relax }}{263}} +\newlabel{ch12:tab:01}{{11.1}{263}} \@writefile{lot}{\contentsline {table}{\numberline {11.2}{\ignorespaces Performances of the parallel CG method on a cluster of 24 CPU cores vs. on a cluster of 12 GPUs.\relax }}{264}} \newlabel{ch12:tab:02}{{11.2}{264}} -\@writefile{lot}{\contentsline {table}{\numberline {11.3}{\ignorespaces Performances of the parallel GMRES method on a cluster 24 CPU cores vs. on cluster of 12 GPUs.\relax }}{265}} -\newlabel{ch12:tab:03}{{11.3}{265}} +\@writefile{lot}{\contentsline {table}{\numberline {11.3}{\ignorespaces Performances of the parallel GMRES method on a cluster 24 CPU cores vs. on cluster of 12 GPUs.\relax }}{264}} +\newlabel{ch12:tab:03}{{11.3}{264}} \newlabel{ch12:eq:20}{{11.20}{265}} \@writefile{lof}{\contentsline {figure}{\numberline {11.6}{\ignorespaces Parallel generation of a large sparse matrix by four computing nodes.\relax }}{266}} \newlabel{ch12:fig:06}{{11.6}{266}} -\@writefile{lot}{\contentsline {table}{\numberline {11.4}{\ignorespaces Main characteristics of sparse banded matrices generated from those of the Davis's collection.\relax }}{267}} -\newlabel{ch12:tab:04}{{11.4}{267}} +\@writefile{lot}{\contentsline {table}{\numberline {11.4}{\ignorespaces Main characteristics of sparse banded matrices generated from those of the Davis collection.\relax }}{266}} +\newlabel{ch12:tab:04}{{11.4}{266}} \@writefile{lot}{\contentsline {table}{\numberline {11.5}{\ignorespaces Performances of the parallel CG method for solving linear systems associated to sparse banded matrices on a cluster of 24 CPU cores vs. on a cluster of 12 GPUs.\relax }}{267}} \newlabel{ch12:tab:05}{{11.5}{267}} -\@writefile{toc}{\contentsline {section}{\numberline {11.5}Hypergraph partitioning}{267}} +\@writefile{lot}{\contentsline {table}{\numberline {11.6}{\ignorespaces Performances of the parallel GMRES method for solving linear systems associated to sparse banded matrices on a cluster of 24 CPU cores vs. on a cluster of 12 GPUs.\relax }}{267}} +\newlabel{ch12:tab:06}{{11.6}{267}} +\@writefile{toc}{\contentsline {section}{\numberline {11.5}Conclusion}{267}} \newlabel{ch12:sec:05}{{11.5}{267}} -\@writefile{lot}{\contentsline {table}{\numberline {11.6}{\ignorespaces Performances of the parallel GMRES method for solving linear systems associated to sparse banded matrices on a cluster of 24 CPU cores vs. on a cluster of 12 GPUs.\relax }}{268}} -\newlabel{ch12:tab:06}{{11.6}{268}} -\@writefile{lot}{\contentsline {table}{\numberline {11.7}{\ignorespaces Main characteristics of sparse five-bands matrices generated from those of the Davis's collection.\relax }}{268}} -\newlabel{ch12:tab:07}{{11.7}{268}} -\@writefile{lof}{\contentsline {figure}{\numberline {11.7}{\ignorespaces Parallel generation of a large sparse five-bands matrix by four computing nodes.\relax }}{269}} -\newlabel{ch12:fig:07}{{11.7}{269}} -\@writefile{lot}{\contentsline {table}{\numberline {11.8}{\ignorespaces Performances of parallel CG solver for solving linear systems associated to sparse five-bands matrices on a cluster of 24 CPU cores vs. on a cluster of 12 GPUs\relax }}{269}} -\newlabel{ch12:tab:08}{{11.8}{269}} -\@writefile{lot}{\contentsline {table}{\numberline {11.9}{\ignorespaces Performances of parallel GMRES solver for solving linear systems associated to sparse five-bands matrices on a cluster of 24 CPU cores vs. on a cluster of 12 GPUs\relax }}{270}} -\newlabel{ch12:tab:09}{{11.9}{270}} -\@writefile{lof}{\contentsline {figure}{\numberline {11.8}{\ignorespaces An example of the hypergraph partitioning of a sparse matrix decomposed between three computing nodes.\relax }}{271}} -\newlabel{ch12:fig:08}{{11.8}{271}} -\@writefile{lot}{\contentsline {table}{\numberline {11.10}{\ignorespaces Performances of the parallel CG solver using hypergraph partitioning for solving linear systems associated to sparse five-bands matrices on a cluster of 24 CPU cores vs. on a cluster of 12 GPU.\relax }}{272}} -\newlabel{ch12:tab:10}{{11.10}{272}} -\@writefile{lot}{\contentsline {table}{\numberline {11.11}{\ignorespaces Performances of the parallel GMRES solver using hypergraph partitioning for solving linear systems associated to sparse five-bands matrices on a cluster of 24 CPU cores vs. on a cluster of 12 GPU.\relax }}{273}} -\newlabel{ch12:tab:11}{{11.11}{273}} -\@writefile{lot}{\contentsline {table}{\numberline {11.12}{\ignorespaces The total communication volume between 12 GPU computing nodes without and with the hypergraph partitioning method.\relax }}{274}} -\newlabel{ch12:tab:12}{{11.12}{274}} -\newlabel{ch12:fig:09.01}{{11.9(a)}{275}} -\newlabel{sub@ch12:fig:09.01}{{(a)}{275}} -\newlabel{ch12:fig:09.02}{{11.9(b)}{275}} -\newlabel{sub@ch12:fig:09.02}{{(b)}{275}} -\@writefile{lof}{\contentsline {figure}{\numberline {11.9}{\ignorespaces Weak-scaling of the parallel CG and GMRES solvers on a GPU cluster for solving large sparse linear systems.\relax }}{275}} -\@writefile{lof}{\contentsline {subfigure}{\numberline{(a)}{\ignorespaces {Sparse band matrices}}}{275}} -\@writefile{lof}{\contentsline {subfigure}{\numberline{(b)}{\ignorespaces {Sparse five-bands matrices}}}{275}} -\newlabel{ch12:fig:09}{{11.9}{275}} -\@writefile{toc}{\contentsline {section}{\numberline {11.6}Conclusion}{275}} -\newlabel{ch12:sec:06}{{11.6}{275}} -\@writefile{toc}{\contentsline {section}{Bibliography}{276}} +\@writefile{toc}{\contentsline {section}{Bibliography}{268}} \@setckpt{Chapters/chapter12/ch12}{ -\setcounter{page}{278} -\setcounter{equation}{25} +\setcounter{page}{270} +\setcounter{equation}{22} \setcounter{enumi}{4} \setcounter{enumii}{0} \setcounter{enumiii}{0} -\setcounter{enumiv}{15} +\setcounter{enumiv}{10} \setcounter{footnote}{0} \setcounter{mpfootnote}{0} \setcounter{part}{1} \setcounter{chapter}{11} -\setcounter{section}{6} +\setcounter{section}{5} \setcounter{subsection}{0} \setcounter{subsubsection}{0} \setcounter{paragraph}{0} \setcounter{subparagraph}{0} -\setcounter{figure}{9} -\setcounter{table}{12} +\setcounter{figure}{6} +\setcounter{table}{6} \setcounter{numauthors}{0} \setcounter{parentequation}{46} \setcounter{subfigure}{0} diff --git a/BookGPU/Chapters/chapter12/ch12.tex b/BookGPU/Chapters/chapter12/ch12.tex index eaa4f9c..7495859 100755 --- a/BookGPU/Chapters/chapter12/ch12.tex +++ b/BookGPU/Chapters/chapter12/ch12.tex @@ -38,11 +38,8 @@ traditional CPUs. In Section~\ref{ch12:sec:02}, we describe the general principle of two well-known iterative methods: the conjugate gradient method and the generalized minimal residual method. In Section~\ref{ch12:sec:03}, we give the main key points of the parallel implementation of both methods on a cluster of -GPUs. Then, in Section~\ref{ch12:sec:04}, we present the experimental results obtained on a -CPU cluster and on a GPU cluster, for solving sparse linear systems associated to matrices -of different structures. Finally, in Section~\ref{ch12:sec:05}, we apply the hypergraph partitioning -technique to reduce the total communication volume between the computing nodes and, thus, -to improve the execution times of the parallel algorithms of both iterative methods. +GPUs. Finally, in Section~\ref{ch12:sec:04}, we present the experimental results obtained on a +CPU cluster and on a GPU cluster, for solving large sparse linear systems. %%--------------------------%% @@ -200,7 +197,7 @@ $maxiter$ are reached. %%****************%% \subsection{GMRES method} \label{ch12:sec:02.02} -The iterative GMRES method is developed by Saad and Schultz in 1986~\cite{ch12:ref3} as a generalization +The iterative GMRES method was developed by Saad and Schultz in 1986~\cite{ch12:ref3} as a generalization of the minimum residual method MINRES~\cite{ch12:ref4}\index{Iterative~method!MINRES}. Indeed, GMRES can be applied for solving symmetric or nonsymmetric linear systems. @@ -231,7 +228,7 @@ V_k = \{v_1, v_2,\ldots,v_k\}, & \forall k>1, v_k=A^{k-1}v_1, \end{equation} and \begin{equation} -V_k A = V_{k+1} \bar{H}_k. +A V_k = V_{k+1} \bar{H}_k. \label{ch12:eq:15} \end{equation} @@ -512,12 +509,6 @@ Tesla C1060 GPU contains $240$ cores running at $1.3$GHz and providing a global a memory bandwidth of $102$GB/s. Figure~\ref{ch12:fig:04} shows the general scheme of the GPU cluster\index{GPU~cluster} that we used in the experimental tests. -\begin{figure} -\centerline{\includegraphics[scale=0.25]{Chapters/chapter12/figures/cluster}} -\caption{General scheme of the GPU cluster of tests composed of six machines, each with two GPUs.} -\label{ch12:fig:04} -\end{figure} - Linux cluster version 2.6.39 OS is installed on CPUs. C programming language is used for coding the parallel algorithms of both methods on the GPU cluster. CUDA version 4.0~\cite{ch12:ref9} is used for programming GPUs, using CUBLAS library~\cite{ch12:ref6} to deal with vector operations @@ -525,6 +516,12 @@ in GPUs and, finally, MPI routines of OpenMPI 1.3.3 are used to carry out the co CPU cores. Indeed, the experiments are done on a cluster of $12$ computing nodes, where each node is managed by a MPI process and it is composed of one CPU core and one GPU card. +\begin{figure}[!h] +\centerline{\includegraphics[scale=0.25]{Chapters/chapter12/figures/cluster}} +\caption{General scheme of the GPU cluster of tests composed of six machines, each with two GPUs.} +\label{ch12:fig:04} +\end{figure} + All tests are made on double-precision floating point operations. The parameters of both linear solvers are initialized as follows: the residual tolerance threshold $\varepsilon=10^{-12}$, the maximum number of iterations $maxiter=500$, the right-hand side $b$ is filled with $1.0$ and the @@ -536,17 +533,9 @@ not too ill-conditioned matrices. In the GPU computing, the size of thread block threads. Finally, the performance results, presented hereafter, are obtained from the mean value over $10$ executions of the same parallel linear solver and for the same input data. -To get more realistic results, we tested the CG and GMRES algorithms on sparse matrices of the Davis's -collection~\cite{ch12:ref10}, that arise in a wide spectrum of real-world applications. We chose six -symmetric sparse matrices and six nonsymmetric ones from this collection. In Figure~\ref{ch12:fig:05}, -we show structures of these matrices and in Table~\ref{ch12:tab:01} we present their main characteristics -which are the number of rows, the total number of nonzero values (nnz) and the maximal bandwidth. In -the present chapter, the bandwidth of a sparse matrix is defined as the number of matrix columns separating -the first and the last nonzero value on a matrix row. - \begin{figure} \centerline{\includegraphics[scale=0.30]{Chapters/chapter12/figures/matrices}} -\caption{Sketches of sparse matrices chosen from the Davis's collection.} +\caption{Sketches of sparse matrices chosen from the Davis collection.} \label{ch12:fig:05} \end{figure} @@ -580,11 +569,18 @@ the first and the last nonzero value on a matrix row. & torso3 & $259,156$ & $4,429,042$ & $216,854$ \\ \hline \end{tabular} -\vspace{0.5cm} -\caption{Main characteristics of sparse matrices chosen from the Davis's collection.} +\caption{Main characteristics of sparse matrices chosen from the Davis collection.} \label{ch12:tab:01} \end{table} +To get more realistic results, we tested the CG and GMRES algorithms on sparse matrices of the Davis +collection~\cite{ch12:ref10}, that arise in a wide spectrum of real-world applications. We chose six +symmetric sparse matrices and six nonsymmetric ones from this collection. In Figure~\ref{ch12:fig:05}, +we show structures of these matrices and in Table~\ref{ch12:tab:01} we present their main characteristics +which are the number of rows, the total number of nonzero values (nnz) and the maximal bandwidth. In +the present chapter, the bandwidth of a sparse matrix is defined as the number of matrix columns separating +the first and the last nonzero value on a matrix row. + \begin{table} \begin{center} \begin{tabular}{|c|c|c|c|c|c|c|} @@ -608,7 +604,7 @@ thermal2 & $1.172s$ & $0.622s$ & $1.88$ & $ \end{center} \end{table} -\begin{table}[!h] +\begin{table} \begin{center} \begin{tabular}{|c|c|c|c|c|c|c|} \hline @@ -675,19 +671,19 @@ of the solution $x^{gpu}$. Thus, we can see that the solutions obtained on the G were computed with a sufficient accuracy (about $10^{-10}$) and they are, more or less, equivalent to those computed on the CPU cluster with a small difference ranging from $10^{-10}$ to $10^{-26}$. However, we can notice from the relative gains $\tau$ that is not interesting to use multiple -GPUs for solving small sparse linear systems. in fact, a small sparse matrix does not allow to +GPUs for solving small sparse linear systems. In fact, a small sparse matrix does not allow to maximize utilization of GPU cores. In addition, the communications required to synchronize the computations over the cluster increase the idle times of GPUs and slow down further the parallel computations. Consequently, in order to test the performances of the parallel solvers, we developed in C programming -language a generator of large sparse matrices. This generator takes a matrix from the Davis's collection~\cite{ch12:ref10} +language a generator of large sparse matrices. This generator takes a matrix from the Davis collection~\cite{ch12:ref10} as an initial matrix to construct large sparse matrices exceeding ten million of rows. It must be executed in parallel by the MPI processes of the computing nodes, so that each process could construct its sparse sub-matrix. In first experimental tests, we are focused on sparse matrices having a banded structure, because they are those arise in the most of numerical problems. So to generate the global sparse matrix, each MPI process constructs its sub-matrix by performing several copies of an initial sparse matrix chosen -from the Davis's collection. Then, it puts all these copies on the main diagonal of the global matrix +from the Davis collection. Then, it puts all these copies on the main diagonal of the global matrix (see Figure~\ref{ch12:fig:06}). Moreover, the empty spaces between two successive copies in the main diagonal are filled with sub-copies (left-copy and right-copy in Figure~\ref{ch12:fig:06}) of the same initial matrix. @@ -729,7 +725,7 @@ initial matrix. & torso3 & $433,795,264$ & $328,757$ \\ \hline \end{tabular} \vspace{0.5cm} -\caption{Main characteristics of sparse banded matrices generated from those of the Davis's collection.} +\caption{Main characteristics of sparse banded matrices generated from those of the Davis collection.} \label{ch12:tab:04} \end{table} @@ -746,7 +742,7 @@ CG method is characterized by a better convergence\index{Convergence} rate and a time of an iteration than those of the GMRES method. Moreover, an iteration of the parallel GMRES method requires more data exchanges between computing nodes compared to the parallel CG method. -\begin{table}[!h] +\begin{table} \begin{center} \begin{tabular}{|c|c|c|c|c|c|c|} \hline @@ -770,7 +766,7 @@ on a cluster of 12 GPUs.} \end{center} \end{table} -\begin{table}[!h] +\begin{table} \begin{center} \begin{tabular}{|c|c|c|c|c|c|c|} \hline @@ -806,390 +802,11 @@ on a cluster of 12 GPUs.} \end{center} \end{table} - %%--------------------------%% %% SECTION 5 %% %%--------------------------%% -\section{Hypergraph partitioning} -\label{ch12:sec:05} -In this section, we present the performances of both parallel CG and GMRES solvers for solving linear -systems associated to sparse matrices having large bandwidths. Indeed, we are interested on sparse -matrices having the nonzero values distributed along their bandwidths. - -\begin{figure} -\centerline{\includegraphics[scale=0.22]{Chapters/chapter12/figures/generation_1}} -\caption{Parallel generation of a large sparse five-bands matrix by four computing nodes.} -\label{ch12:fig:07} -\end{figure} - -\begin{table}[!h] -\begin{center} -\begin{tabular}{|c|c|c|c|} -\hline -{\bf Matrix type} & {\bf Matrix name} & {\bf \# nnz} & {\bf 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}{*}{Nonsymmetric} & 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 sparse five-bands matrices generated from those of the Davis's collection.} -\label{ch12:tab:07} -\end{center} -\end{table} - -We have developed in C programming language a generator of large sparse matrices -having five bands distributed along their bandwidths (see Figure~\ref{ch12:fig:07}). -The principle of this generator is equivalent to that in Section~\ref{ch12:sec:04}. -However, the copies performed on the initial matrix (chosen from the Davis's collection) -are placed on the main diagonal and on four off-diagonals, two on the right and two -on the left of the main diagonal. Figure~\ref{ch12:fig:07} shows an example of a -generation of a sparse five-bands matrix by four computing nodes. Table~\ref{ch12:tab:07} -shows the main characteristics of sparse five-bands matrices generated from those -presented in Table~\ref{ch12:tab:01} and associated to linear systems of $25$ million -unknown values. - -\begin{table}[!h] -\begin{center} -\begin{tabular}{|c|c|c|c|c|c|c|} -\hline -{\bf Matrix} & $\mathbf{Time_{cpu}}$ & $\mathbf{Time_{gpu}}$ & $\mathbf{\tau}$ & $\mathbf{\# iter.}$ & $\mathbf{prec.}$ & $\mathbf{\Delta}$ \\ \hline \hline - -2cubes\_sphere & $6.041s$ & $3.338s$ & $1.81$ & $30$ & $6.77e$-$11$ & $3.25e$-$19$ \\ - -ecology2 & $1.404s$ & $1.301s$ & $1.08$ & $13$ & $5.22e$-$11$ & $2.17e$-$18$ \\ - -finan512 & $1.822s$ & $1.299s$ & $1.40$ & $12$ & $3.52e$-$11$ & $3.47e$-$18$ \\ - -G3\_circuit & $2.331s$ & $2.129s$ & $1.09$ & $15$ & $1.36e$-$11$ & $5.20e$-$18$ \\ - -shallow\_water2 & $0.541s$ & $0.504s$ & $1.07$ & $6$ & $2.12e$-$16$ & $5.05e$-$28$ \\ - -thermal2 & $2.549s$ & $1.705s$ & $1.49$ & $14$ & $2.36e$-$10$ & $5.20e$-$18$ \\ \hline -\end{tabular} -\caption{Performances of parallel CG solver for solving linear systems associated to sparse five-bands matrices -on a cluster of 24 CPU cores vs. on a cluster of 12 GPUs} -\label{ch12:tab:08} -\end{center} -\end{table} - -\begin{table} -\begin{center} -\begin{tabular}{|c|c|c|c|c|c|c|} -\hline -{\bf Matrix} & $\mathbf{Time_{cpu}}$ & $\mathbf{Time_{gpu}}$ & $\mathbf{\tau}$ & $\mathbf{\# iter.}$ & $\mathbf{prec.}$ & $\mathbf{\Delta}$ \\ \hline \hline - -2cubes\_sphere & $15.963s$ & $7.250s$ & $2.20$ & $58$ & $6.23e$-$16$ & $3.25e$-$19$ \\ - -ecology2 & $3.549s$ & $2.176s$ & $1.63$ & $21$ & $4.78e$-$15$ & $1.06e$-$15$ \\ - -finan512 & $3.862s$ & $1.934s$ & $1.99$ & $17$ & $3.21e$-$14$ & $8.43e$-$17$ \\ - -G3\_circuit & $4.636s$ & $2.811s$ & $1.65$ & $22$ & $1.08e$-$14$ & $1.77e$-$16$ \\ - -shallow\_water2 & $2.738s$ & $1.539s$ & $1.78$ & $17$ & $5.54e$-$23$ & $3.82e$-$26$ \\ - -thermal2 & $5.017s$ & $2.587s$ & $1.94$ & $21$ & $8.25e$-$14$ & $4.34e$-$18$ \\ \hline \hline - -cage13 & $9.315s$ & $3.227s$ & $2.89$ & $26$ & $3.38e$-$13$ & $2.08e$-$16$ \\ - -crashbasis & $35.980s$ & $14.770s$ & $2.43$ & $127$ & $1.17e$-$12$ & $1.56e$-$17$ \\ - -FEM\_3D\_thermal2 & $24.611s$ & $7.749s$ & $3.17$ & $64$ & $3.87e$-$11$ & $2.84e$-$14$ \\ - -language & $16.859s$ & $9.697s$ & $1.74$ & $89$ & $2.17e$-$12$ & $1.70e$-$12$ \\ - -poli\_large & $10.200s$ & $6.534s$ & $1.56$ & $69$ & $5.14e$-$13$ & $1.63e$-$13$ \\ - -torso3 & $49.074s$ & $19.397s$ & $2.53$ & $175$ & $2.69e$-$12$ & $2.77e$-$16$ \\ \hline -\end{tabular} -\caption{Performances of parallel GMRES solver for solving linear systems associated to sparse five-bands matrices -on a cluster of 24 CPU cores vs. on a cluster of 12 GPUs} -\label{ch12:tab:09} -\end{center} -\end{table} - -Tables~\ref{ch12:tab:08} and~\ref{ch12:tab:09} shows the performances of the parallel -CG and GMRES solvers, respectively, obtained on a cluster of $24$ CPU cores and on a -cluster of $12$ GPUs. The linear systems solved in these tables are associated to the -sparse five-bands matrices presented on Table~\ref{ch12:tab:07}. We can notice from -both Tables~\ref{ch12:tab:08} and~\ref{ch12:tab:09} that using a GPU cluster is not -efficient for solving these kind of sparse linear systems\index{Sparse~linear~system}. -We can see that the execution times obtained on the GPU cluster are almost equivalent -to those obtained on the CPU cluster (see the relative gains presented in column~$4$ -of each table). This is due to the large number of communications necessary to synchronize -the computations over the cluster. Indeed, the naive partitioning, row-by-row or column-by-column, -of sparse matrices having large bandwidths can link a computing node to many neighbors -and then generate a large number of data dependencies between these computing nodes in -the cluster. - -Therefore, we have chosen to use a hypergraph partitioning method\index{Hypergraph}, -which is well-suited to numerous kinds of sparse matrices~\cite{ch12:ref11}. Indeed, -it can well model the communications between the computing nodes, particularly in the -case of nonsymmetric and irregular matrices, and it gives good reduction of the total -communication volume. In contrast, it is an expensive operation in terms of execution -time and memory space. - -The sparse matrix $A$ of the linear system to be solved is modeled as a hypergraph -$\mathcal{H}=(\mathcal{V},\mathcal{E})$\index{Hypergraph} as follows: -\begin{itemize} -\item each matrix row $\{i\}_{0\leq i0$ an approximate solution $x_k$ -which, gradually, converges to the exact solution $x^{*}$ as follows: +where $A\in\mathbb{R}^{n\times n}$ is a sparse nonsingular square matrix, $x\in\mathbb{R}^{n}$ +is the solution vector, $b\in\mathbb{R}^{n}$ is the right-hand side and $n\in\mathbb{N}$ is a +large integer number. + +The iterative methods\index{Iterative~method} for solving the large sparse linear system~(\ref{ch12:eq:01}) +proceed by successive iterations of a same block of elementary operations, during which an +infinite number of approximate solutions $\{x_k\}_{k\geq 0}$ are computed. Indeed, from an +initial guess $x_0$, an iterative method determines at each iteration $k>0$ an approximate +solution $x_k$ which, gradually, converges to the exact solution $x^{*}$ as follows: \begin{equation} x^{*}=\lim\limits_{k\to\infty}x_{k}=A^{-1}b. -\label{eq:02} +\label{ch12:eq:02} \end{equation} -The number of iterations necessary to reach the exact solution $x^{*}$ is not known beforehand and can be infinite. In -practice, an iterative method often finds an approximate solution $\tilde{x}$ after a fixed number of iterations and/or -when a given convergence criterion is satisfied as follows: +The number of iterations necessary to reach the exact solution $x^{*}$ is not known beforehand +and can be infinite. In practice, an iterative method often finds an approximate solution $\tilde{x}$ +after a fixed number of iterations and/or when a given convergence criterion\index{Convergence} +is satisfied as follows: \begin{equation} \|b-A\tilde{x}\| < \varepsilon, -\label{eq:03} +\label{ch12:eq:03} \end{equation} -where $\varepsilon<1$ is the required convergence tolerance threshold. - -Some of the most iterative methods that have proven their efficiency for solving large sparse linear systems are those -called \textit{Krylov sub-space methods}~\cite{ref1}. In the present chapter, we describe two Krylov methods which are -widely used: the conjugate gradient method (CG) and the generalized minimal residual method (GMRES). In practice, the -Krylov sub-space methods are usually used with preconditioners that allow to improve their convergence. So, in what -follows, the CG and GMRES methods are used for solving the left-preconditioned sparse linear system: +where $\varepsilon<1$ is the required convergence tolerance threshold\index{Convergence!Tolerance~threshold}. + +Some of the most iterative methods that have proven their efficiency for solving large sparse +linear systems are those called \textit{Krylov subspace methods}~\cite{ch12:ref1}\index{Iterative~method!Krylov~subspace}. +In the present chapter, we describe two Krylov methods which are widely used: the conjugate +gradient method (CG) and the generalized minimal residual method (GMRES). In practice, the +Krylov subspace methods are usually used with preconditioners that allow to improve their +convergence. So, in what follows, the CG and GMRES methods are used for solving the left-preconditioned\index{Sparse~linear~system!Preconditioned} +sparse linear system: \begin{equation} M^{-1}Ax=M^{-1}b, -\label{eq:11} +\label{ch12:eq:11} \end{equation} where $M$ is the preconditioning matrix. + %%****************%% %%****************%% \subsection{CG method} -\label{sec:02.01} -The conjugate gradient method is initially developed by Hestenes and Stiefel in 1952~\cite{ref2}. It is one of the well -known iterative method for solving large sparse linear systems. In addition, it can be adapted for solving nonlinear -equations and optimization problems. However, it can only be applied to problems with positive definite symmetric matrices. - -The main idea of the CG method is the computation of a sequence of approximate solutions $\{x_k\}_{k\geq 0}$ in a Krylov -sub-space of order $k$ as follows: +\label{ch12:sec:02.01} +The conjugate gradient method is initially developed by Hestenes and Stiefel in 1952~\cite{ch12:ref2}. +It is one of the well known iterative method for solving large sparse linear systems. In addition, it +can be adapted for solving nonlinear equations and optimization problems. However, it can only be applied +to problems with positive definite symmetric matrices. + +The main idea of the CG method\index{Iterative~method!CG} is the computation of a sequence of approximate +solutions $\{x_k\}_{k\geq 0}$ in a Krylov subspace\index{Iterative~method!Krylov~subspace} of order $k$ as +follows: \begin{equation} x_k \in x_0 + \mathcal{K}_k(A,r_0), -\label{eq:04} +\label{ch12:eq:04} \end{equation} -such that the Galerkin condition must be satisfied: +such that the Galerkin condition\index{Galerkin~condition} must be satisfied: \begin{equation} r_k \bot \mathcal{K}_k(A,r_0), -\label{eq:05} +\label{ch12:eq:05} \end{equation} -where $x_0$ is the initial guess, $r_k=b-Ax_k$ is the residual of the computed solution $x_k$ and $\mathcal{K}_k$ the Krylov -sub-space of order $k$: \[\mathcal{K}_k(A,r_0) \equiv\text{span}\{r_0, Ar_0, A^2r_0,\ldots, A^{k-1}r_0\}.\] +where $x_0$ is the initial guess, $r_k=b-Ax_k$ is the residual of the computed solution $x_k$ and $\mathcal{K}_k$ +the Krylov subspace of order $k$: \[\mathcal{K}_k(A,r_0) \equiv\text{span}\{r_0, Ar_0, A^2r_0,\ldots, A^{k-1}r_0\}.\] In fact, CG is based on the construction of a sequence $\{p_k\}_{k\in\mathbb{N}}$ of direction vectors in $\mathcal{K}_k$ which are pairwise $A$-conjugate ($A$-orthogonal): \begin{equation} \begin{array}{ll} p_i^T A p_j = 0, & i\neq j. \end{array} -\label{eq:06} +\label{ch12:eq:06} \end{equation} At each iteration $k$, an approximate solution $x_k$ is computed by recurrence as follows: \begin{equation} \begin{array}{ll} x_k = x_{k-1} + \alpha_k p_k, & \alpha_k\in\mathbb{R}. \end{array} -\label{eq:07} +\label{ch12:eq:07} \end{equation} Consequently, the residuals $r_k$ are computed in the same way: \begin{equation} r_k = r_{k-1} - \alpha_k A p_k. -\label{eq:08} +\label{ch12:eq:08} \end{equation} -In the case where all residuals are nonzero, the direction vectors $p_k$ can be determined so that the following recurrence -holds: +In the case where all residuals are nonzero, the direction vectors $p_k$ can be determined so that +the following recurrence holds: \begin{equation} \begin{array}{lll} p_0=r_0, & p_k=r_k+\beta_k p_{k-1}, & \beta_k\in\mathbb{R}. \end{array} -\label{eq:09} +\label{ch12:eq:09} \end{equation} -Moreover, the scalars $\{\alpha_k\}_{k>0}$ are chosen so as to minimize the $A$-norm error $\|x^{*}-x_k\|_A$ over the Krylov -sub-space $\mathcal{K}_{k}$ and the scalars $\{\beta_k\}_{k>0}$ are chosen so as to ensure that the direction vectors are -pairwise $A$-conjugate. So, the assumption that matrix $A$ is symmetric and the recurrences~(\ref{eq:08}) and~(\ref{eq:09}) -allow to deduce that: +Moreover, the scalars $\{\alpha_k\}_{k>0}$ are chosen so as to minimize the $A$-norm error $\|x^{*}-x_k\|_A$ +over the Krylov subspace $\mathcal{K}_{k}$ and the scalars $\{\beta_k\}_{k>0}$ are chosen so as to ensure +that the direction vectors are pairwise $A$-conjugate. So, the assumption that matrix $A$ is symmetric and +the recurrences~(\ref{ch12:eq:08}) and~(\ref{ch12:eq:09}) allow to deduce that: \begin{equation} \begin{array}{ll} \alpha_{k}=\frac{r^{T}_{k-1}r_{k-1}}{p_{k}^{T}Ap_{k}}, & \beta_{k}=\frac{r_{k}^{T}r_{k}}{r_{k-1}^{T}r_{k-1}}. \end{array} -\label{eq:10} +\label{ch12:eq:10} \end{equation} \begin{algorithm}[!t] - \SetLine - \linesnumbered Choose an initial guess $x_0$\; $r_{0} = b - A x_{0}$\; $convergence$ = false\; @@ -160,62 +180,70 @@ allow to deduce that: } } \caption{Left-preconditioned CG method} -\label{alg:01} +\label{ch12:alg:01} \end{algorithm} -Algorithm~\ref{alg:01} shows the main key points of the preconditioned CG method. It allows to solve the left-preconditioned -sparse linear system~(\ref{eq:11}). In this algorithm, $\varepsilon$ is the convergence tolerance threshold, $maxiter$ is the maximum -number of iterations and $(\cdot,\cdot)$ defines the dot product between two vectors in $\mathbb{R}^{n}$. At every iteration, a direction -vector $p_k$ is determined, so that it is orthogonal to the preconditioned residual $z_k$ and to the direction vectors $\{p_i\}_{i0}$ in a Krylov sub-space $\mathcal{K}_k$ as follows: +\label{ch12:sec:02.02} +The iterative GMRES method is developed by Saad and Schultz in 1986~\cite{ch12:ref3} as a generalization +of the minimum residual method MINRES~\cite{ch12:ref4}\index{Iterative~method!MINRES}. Indeed, GMRES can +be applied for solving symmetric or nonsymmetric linear systems. + +The main principle of the GMRES method\index{Iterative~method!GMRES} is to find an approximation minimizing +at best the residual norm. In fact, GMRES computes a sequence of approximate solutions $\{x_k\}_{k>0}$ in +a Krylov subspace\index{Iterative~method!Krylov~subspace} $\mathcal{K}_k$ as follows: \begin{equation} \begin{array}{ll} x_k \in x_0 + \mathcal{K}_k(A, v_1),& v_1=\frac{r_0}{\|r_0\|_2}, \end{array} -\label{eq:12} +\label{ch12:eq:12} \end{equation} -so that the Petrov-Galerkin condition is satisfied: +so that the Petrov-Galerkin condition\index{Petrov-Galerkin~condition} is satisfied: \begin{equation} \begin{array}{ll} r_k \bot A \mathcal{K}_k(A, v_1). \end{array} -\label{eq:13} +\label{ch12:eq:13} \end{equation} -GMRES uses the Arnoldi process~\cite{ref5} to construct an orthonormal basis $V_k$ for the Krylov sub-space $\mathcal{K}_k$ -and an upper Hessenberg matrix $\bar{H}_k$ of order $(k+1)\times k$: +GMRES uses the Arnoldi process~\cite{ch12:ref5}\index{Iterative~method!Arnoldi~process} to construct an +orthonormal basis $V_k$ for the Krylov subspace $\mathcal{K}_k$ and an upper Hessenberg matrix\index{Hessenberg~matrix} +$\bar{H}_k$ of order $(k+1)\times k$: \begin{equation} \begin{array}{ll} V_k = \{v_1, v_2,\ldots,v_k\}, & \forall k>1, v_k=A^{k-1}v_1, \end{array} -\label{eq:14} +\label{ch12:eq:14} \end{equation} and \begin{equation} V_k A = V_{k+1} \bar{H}_k. -\label{eq:15} +\label{ch12:eq:15} \end{equation} -Then, at each iteration $k$, an approximate solution $x_k$ is computed in the Krylov sub-space $\mathcal{K}_k$ spanned by $V_k$ -as follows: +Then, at each iteration $k$, an approximate solution $x_k$ is computed in the Krylov subspace $\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:16} +\label{ch12:eq:16} \end{equation} -From both formulas~(\ref{eq:15}) and~(\ref{eq:16}) and $r_k=b-Ax_k$, we can deduce that: +From both formulas~(\ref{ch12:eq:15}) and~(\ref{ch12:eq:16}) and $r_k=b-Ax_k$, we can deduce that: \begin{equation} \begin{array}{lll} r_{k} & = & b - A (x_{0} + V_{k}y) \\ @@ -223,34 +251,34 @@ From both formulas~(\ref{eq:15}) and~(\ref{eq:16}) and $r_k=b-Ax_k$, we can dedu & = & \beta v_{1} - V_{k+1}\bar{H}_{k}y \\ & = & V_{k+1}(\beta e_{1} - \bar{H}_{k}y), \end{array} -\label{eq:17} +\label{ch12:eq:17} \end{equation} -such that $\beta=\|r_0\|_2$ and $e_1=(1,0,\cdots,0)$ is the first vector of the canonical basis of $\mathbb{R}^k$. So, -the vector $y$ is chosen in $\mathbb{R}^k$ so as to minimize at best the Euclidean norm of the residual $r_k$. Consequently, -a linear least-squares problem of size $k$ is solved: +such that $\beta=\|r_0\|_2$ and $e_1=(1,0,\cdots,0)$ is the first vector of the canonical basis of +$\mathbb{R}^k$. So, the vector $y$ is chosen in $\mathbb{R}^k$ so as to minimize at best the Euclidean +norm of the residual $r_k$. Consequently, a linear least-squares problem of size $k$ is 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}. -\label{eq:18} +\label{ch12:eq:18} \end{equation} -The QR factorization of matrix $\bar{H}_k$ is used to compute the solution of this problem by using Givens rotations~\cite{ref1,ref3}, -such that: +The QR factorization of matrix $\bar{H}_k$ is used to compute the solution of this problem by using +Givens rotations~\cite{ch12:ref1,ch12:ref3}, 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} -\label{eq:19} +\label{ch12:eq:19} \end{equation} where $Q_kQ_k^T=I_k$ and $R_k$ is an upper triangular matrix. -The GMRES method computes an approximate solution with a sufficient precision after, at most, $n$ iterations ($n$ is the size of the -sparse linear system to be solved). However, the GMRES algorithm must construct and store in the memory an orthonormal basis $V_k$ whose -size is proportional to the number of iterations required to achieve the convergence. Then, to avoid a huge memory storage, the GMRES -method must be restarted at each $m$ iterations, such that $m$ is very small ($m\ll n$), and with $x_m$ as the initial guess to the -next iteration. This allows to limit the size of the basis $V$ to $m$ orthogonal vectors. +The GMRES method computes an approximate solution with a sufficient precision after, at most, $n$ +iterations ($n$ is the size of the sparse linear system to be solved). However, the GMRES algorithm +must construct and store in the memory an orthonormal basis $V_k$ whose size is proportional to the +number of iterations required to achieve the convergence. Then, to avoid a huge memory storage, the +GMRES method must be restarted at each $m$ iterations, such that $m$ is very small ($m\ll n$), and +with $x_m$ as the initial guess to the next iteration. This allows to limit the size of the basis +$V$ to $m$ orthogonal vectors. \begin{algorithm}[!t] - \SetLine - \linesnumbered Choose an initial guess $x_0$\; $convergence$ = false\; $k = 1$\; @@ -281,190 +309,245 @@ next iteration. This allows to limit the size of the basis $V$ to $m$ orthogonal } } \caption{Left-preconditioned GMRES method with restarts} -\label{alg:02} +\label{ch12:alg:02} \end{algorithm} -Algorithm~\ref{alg:02} shows the main key points of the GMRES method with restarts. It solves the left-preconditioned sparse linear -system~(\ref{eq:11}), such that $M$ is the preconditioning matrix. At each iteration $k$, GMRES uses the Arnoldi process (defined -from line~$7$ to line~$17$) to construct a basis $V_m$ of $m$ orthogonal vectors and an upper Hessenberg matrix $\bar{H}_m$ of size -$(m+1)\times m$. Then, it solves the linear least-squares problem of size $m$ to find the vector $y\in\mathbb{R}^{m}$ which minimizes -at best the residual norm (line~$18$). Finally, it computes an approximate solution $x_m$ in the Krylov sub-space spanned by $V_m$ -(line~$19$). The GMRES algorithm is stopped when the residual norm is sufficiently small ($\|r_m\|_2<\varepsilon$) and/or the maximum -number of iterations ($maxiter$) is reached. +Algorithm~\ref{ch12:alg:02} shows the main key points of the GMRES method with restarts. +It solves the left-preconditioned\index{Sparse~linear~system!Preconditioned} sparse linear +system~(\ref{ch12:eq:11}), such that $M$ is the preconditioning matrix. At each iteration +$k$, GMRES uses the Arnoldi process\index{Iterative~method!Arnoldi~process} (defined from +line~$7$ to line~$17$) to construct a basis $V_m$ of $m$ orthogonal vectors and an upper +Hessenberg matrix\index{Hessenberg~matrix} $\bar{H}_m$ of size $(m+1)\times m$. Then, it +solves the linear least-squares problem of size $m$ to find the vector $y\in\mathbb{R}^{m}$ +which minimizes at best the residual norm (line~$18$). Finally, it computes an approximate +solution $x_m$ in the Krylov subspace spanned by $V_m$ (line~$19$). The GMRES algorithm is +stopped when the residual norm is sufficiently small ($\|r_m\|_2<\varepsilon$) and/or the +maximum number of iterations\index{Convergence!Maximum~number~of~iterations} ($maxiter$) +is reached. + %%--------------------------%% %% SECTION 3 %% %%--------------------------%% \section{Parallel implementation on a GPU cluster} -\label{sec:03} -In this section, we present the parallel algorithms of both iterative CG and GMRES methods for GPU clusters. -The implementation is performed on a GPU cluster composed of different computing nodes, such that each node -is a CPU core managed by a MPI process and equipped with a GPU card. The parallelization of these algorithms -is carried out by using the MPI communication routines between the GPU computing nodes and the CUDA programming -environment inside each node. In what follows, the algorithms of the iterative methods are called iterative -solvers. +\label{ch12:sec:03} +In this section, we present the parallel algorithms of both iterative CG\index{Iterative~method!CG} +and GMRES\index{Iterative~method!GMRES} methods for GPU clusters. The implementation is performed on +a GPU cluster composed of different computing nodes, such that each node is a CPU core managed by a +MPI process and equipped with a GPU card. The parallelization of these algorithms is carried out by +using the MPI communication routines between the GPU computing nodes\index{Computing~node} and the +CUDA programming environment inside each node. In what follows, the algorithms of the iterative methods +are called iterative solvers. + %%****************%% %%****************%% \subsection{Data partitioning} -\label{sec:03.01} -The parallel solving of the large sparse linear system~(\ref{eq:11}) requires a data partitioning between the computing -nodes of the GPU cluster. Let $p$ denotes the number of the computing nodes on the GPU cluster. The partitioning operation -consists in the decomposition of the vectors and matrices, involved in the iterative solver, in $p$ portions. Indeed, this -operation allows to assign to each computing node $i$: -\begin{itemize*} +\label{ch12:sec:03.01} +The parallel solving of the large sparse linear system~(\ref{ch12:eq:11}) requires a data partitioning +between the computing nodes of the GPU cluster. Let $p$ denotes the number of the computing nodes on the +GPU cluster. The partitioning operation consists in the decomposition of the vectors and matrices, involved +in the iterative solver, in $p$ portions. Indeed, this operation allows to assign to each computing node +$i$: +\begin{itemize} \item a portion of size $\frac{n}{p}$ elements of each vector, \item a sparse rectangular sub-matrix $A_i$ of size $(\frac{n}{p},n)$ and, \item a square preconditioning sub-matrix $M_i$ of size $(\frac{n}{p},\frac{n}{p})$, -\end{itemize*} -where $n$ is the size of the sparse linear system to be solved. In the first instance, we perform a naive row-wise partitioning -(decomposition row-by-row) on the data of the sparse linear systems to be solved. Figure~\ref{fig:01} shows an example of a row-wise -data partitioning between four computing nodes of a sparse linear system (sparse matrix $A$, solution vector $x$ and right-hand -side $b$) of size $16$ unknown values. +\end{itemize} +where $n$ is the size of the sparse linear system to be solved. In the first instance, we perform a naive +row-wise partitioning (decomposition row-by-row) on the data of the sparse linear systems to be solved. +Figure~\ref{ch12:fig:01} shows an example of a row-wise data partitioning between four computing nodes +of a sparse linear system (sparse matrix $A$, solution vector $x$ and right-hand side $b$) of size $16$ +unknown values. \begin{figure} \centerline{\includegraphics[scale=0.35]{Chapters/chapter12/figures/partition}} \caption{A data partitioning of the sparse matrix $A$, the solution vector $x$ and the right-hand side $b$ into four portions.} -\label{fig:01} +\label{ch12:fig:01} \end{figure} + %%****************%% %%****************%% \subsection{GPU computing} -\label{sec:03.02} -After the partitioning operation, all the data involved from this operation must be transferred from the CPU memories to the GPU -memories, in order to be processed by GPUs. We use two functions of the CUBLAS library (CUDA Basic Linear Algebra Subroutines), -developed by Nvidia~\cite{ref6}: \verb+cublasAlloc()+ for the memory allocations on GPUs and \verb+cublasSetVector()+ for the -memory copies from the CPUs to the GPUs. - -An efficient implementation of CG and GMRES solvers on a GPU cluster requires to determine all parts of their codes that can be -executed in parallel and, thus, take advantage of the GPU acceleration. As many Krylov sub-space methods, the CG and GMRES methods -are mainly based on arithmetic operations dealing with vectors or matrices: sparse matrix-vector multiplications, scalar-vector -multiplications, dot products, Euclidean norms, AXPY operations ($y\leftarrow ax+y$ where $x$ and $y$ are vectors and $a$ is a -scalar) and so on. These vector operations are often easy to parallelize and they are more efficient on parallel computers when -they work on large vectors. Therefore, all the vector operations used in CG and GMRES solvers must be executed by the GPUs as kernels. - -We use the kernels of the CUBLAS library to compute some vector operations of CG and GMRES solvers. The following kernels of CUBLAS -(dealing with double floating point) are used: \verb+cublasDdot()+ for the dot products, \verb+cublasDnrm2()+ for the Euclidean -norms and \verb+cublasDaxpy()+ for the AXPY operations. For the rest of the data-parallel operations, we code their kernels in CUDA. -In the CG solver, we develop a kernel for the XPAY operation ($y\leftarrow x+ay$) used at line~$12$ in Algorithm~\ref{alg:01}. In the -GMRES solver, we program a kernel for the scalar-vector multiplication (lines~$7$ and~$15$ in Algorithm~\ref{alg:02}), a kernel for -solving the least-squares problem and a kernel for the elements updates of the solution vector $x$. - -The least-squares problem in the GMRES method is solved by performing a QR factorization on the Hessenberg matrix $\bar{H}_m$ with -plane rotations and, then, solving the triangular system by backward substitutions to compute $y$. Consequently, solving the least-squares -problem on the GPU is not interesting. Indeed, the triangular solves are not easy to parallelize and inefficient on GPUs. However, -the least-squares problem to solve in the GMRES method with restarts has, generally, a very small size $m$. Therefore, we develop -an inexpensive kernel which must be executed in sequential by a single CUDA thread. - -The most important operation in CG and GMRES methods is the sparse matrix-vector multiplication (SpMV), because it is often an -expensive operation in terms of execution time and memory space. Moreover, it requires to take care of the storage format of the -sparse matrix in the memory. Indeed, the naive storage, row-by-row or column-by-column, of a sparse matrix can cause a significant -waste of memory space and execution time. In addition, the sparsity nature of the matrix often leads to irregular memory accesses -to read the matrix nonzero values. So, the computation of the SpMV multiplication on GPUs can involve non coalesced accesses to -the global memory, which slows down even more its performances. One of the most efficient compressed storage formats of sparse -matrices on GPUs is HYB format~\cite{ref7}. It is a combination of ELLpack (ELL) and Coordinate (COO) formats. Indeed, 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 accesses and the flexibility of COO which is insensitive to the matrix -structure. Consequently, we use the HYB kernel~\cite{ref8} developed by Nvidia to implement the SpMV multiplication of CG and -GMRES methods on GPUs. Moreover, to avoid the non coalesced accesses to the high-latency global memory, we fill the elements of -the iterate vector $x$ in the cached texture memory. +\label{ch12:sec:03.02} +After the partitioning operation, all the data involved from this operation must be +transferred from the CPU memories to the GPU memories, in order to be processed by +GPUs. We use two functions of the CUBLAS\index{CUBLAS} library (CUDA Basic Linear +Algebra Subroutines), developed by Nvidia~\cite{ch12:ref6}: \verb+cublasAlloc()+ +for the memory allocations on GPUs and \verb+cublasSetVector()+ for the memory +copies from the CPUs to the GPUs. + +An efficient implementation of CG and GMRES solvers on a GPU cluster requires to +determine all parts of their codes that can be executed in parallel and, thus, take +advantage of the GPU acceleration. As many Krylov subspace methods, the CG and GMRES +methods are mainly based on arithmetic operations dealing with vectors or matrices: +sparse matrix-vector multiplications, scalar-vector multiplications, dot products, +Euclidean norms, AXPY operations ($y\leftarrow ax+y$ where $x$ and $y$ are vectors +and $a$ is a scalar) and so on. These vector operations are often easy to parallelize +and they are more efficient on parallel computers when they work on large vectors. +Therefore, all the vector operations used in CG and GMRES solvers must be executed +by the GPUs as kernels. + +We use the kernels of the CUBLAS library to compute some vector operations of CG and +GMRES solvers. The following kernels of CUBLAS (dealing with double floating point) +are used: \verb+cublasDdot()+ for the dot products, \verb+cublasDnrm2()+ for the +Euclidean norms and \verb+cublasDaxpy()+ for the AXPY operations. For the rest of +the data-parallel operations, we code their kernels in CUDA. In the CG solver, we +develop a kernel for the XPAY operation ($y\leftarrow x+ay$) used at line~$12$ in +Algorithm~\ref{ch12:alg:01}. In the GMRES solver, we program a kernel for the scalar-vector +multiplication (lines~$7$ and~$15$ in Algorithm~\ref{ch12:alg:02}), a kernel for +solving the least-squares problem and a kernel for the elements updates of the solution +vector $x$. + +The least-squares problem in the GMRES method is solved by performing a QR factorization +on the Hessenberg matrix\index{Hessenberg~matrix} $\bar{H}_m$ with plane rotations and, +then, solving the triangular system by backward substitutions to compute $y$. Consequently, +solving the least-squares problem on the GPU is not interesting. Indeed, the triangular +solves are not easy to parallelize and inefficient on GPUs. However, the least-squares +problem to solve in the GMRES method with restarts has, generally, a very small size $m$. +Therefore, we develop an inexpensive kernel which must be executed in sequential by a +single CUDA thread. + +The most important operation in CG\index{Iterative~method!CG} and GMRES\index{Iterative~method!GMRES} +methods is the sparse matrix-vector multiplication (SpMV)\index{SpMV~multiplication}, +because it is often an expensive operation in terms of execution time and memory space. +Moreover, it requires to take care of the storage format of the sparse matrix in the +memory. Indeed, the naive storage, row-by-row or column-by-column, of a sparse matrix +can cause a significant waste of memory space and execution time. In addition, the sparsity +nature of the matrix often leads to irregular memory accesses to read the matrix nonzero +values. So, the computation of the SpMV multiplication on GPUs can involve non coalesced +accesses to the global memory, which slows down even more its performances. One of the +most efficient compressed storage formats\index{Compressed~storage~format} of sparse +matrices on GPUs is HYB\index{Compressed~storage~format!HYB} format~\cite{ch12:ref7}. +It is a combination of ELLpack (ELL) and Coordinate (COO) formats. Indeed, it stores +a typical number of nonzero values per row in ELL\index{Compressed~storage~format!ELL} +format and remaining entries of exceptional rows in COO format. It combines the efficiency +of ELL due to the regularity of its memory accesses and the flexibility of COO\index{Compressed~storage~format!COO} +which is insensitive to the matrix structure. Consequently, we use the HYB kernel~\cite{ch12:ref8} +developed by Nvidia to implement the SpMV multiplication of CG and GMRES methods on GPUs. +Moreover, to avoid the non coalesced accesses to the high-latency global memory, we fill +the elements of the iterate vector $x$ in the cached texture memory. + %%****************%% %%****************%% \subsection{Data communications} -\label{sec:03.03} -All the computing nodes of the GPU cluster execute in parallel the same iterative solver (Algorithm~\ref{alg:01} or Algorithm~\ref{alg:02}) -adapted to GPUs, but on their own portions of the sparse linear system: $M^{-1}_iA_ix_i=M^{-1}_ib_i$, $0\leq i