new

[book_gpu.git] / BookGPU / Chapters / chapter12 / ch12.tex

diff --git a/BookGPU/Chapters/chapter12/ch12.tex b/BookGPU/Chapters/chapter12/ch12.tex
index 749585952548f74f180bca570c0db36d0e6f0a63..254c0cbd83c36dbe0cb93101bc1894637ba72bd9 100755 (executable)
--- a/BookGPU/Chapters/chapter12/ch12.tex
+++ b/BookGPU/Chapters/chapter12/ch12.tex
@@ -17,10 +17,10 @@
 %%--------------------------%%
 \section{Introduction}
 \label{ch12:sec:01}
 %%--------------------------%%
 \section{Introduction}
 \label{ch12:sec:01}

-The sparse linear systems are used to model many scientific and industrial problems,
+Sparse linear systems are used to model many scientific and industrial problems,

 such as the environmental simulations or the industrial processing of the complex or
 non-Newtonian fluids. Moreover, the resolution of these problems often involves the
 such as the environmental simulations or the industrial processing of the complex or
 non-Newtonian fluids. Moreover, the resolution of these problems often involves the

-solving of such linear systems which is considered as the most expensive process in
+solving of such linear systems which are considered as the most expensive process in

 terms of execution time and memory space. Therefore, solving sparse linear systems
 must be as efficient as possible in order to deal with problems of ever increasing
 size.
 terms of execution time and memory space. Therefore, solving sparse linear systems
 must be as efficient as possible in order to deal with problems of ever increasing
 size.


@@ -28,10 +28,10 @@ size.
 There are, in the jargon of numerical analysis, different methods of solving sparse
 linear systems that can be classified in two classes: the direct and iterative methods.
 However, the iterative methods are often more suitable than their counterpart, direct
 There are, in the jargon of numerical analysis, different methods of solving sparse
 linear systems that can be classified in two classes: the direct and iterative methods.
 However, the iterative methods are often more suitable than their counterpart, direct

-methods, for solving these systems. Indeed, they are less memory consuming and easier
+methods, to  solve these systems. Indeed, they are less memory consuming and easier

 to parallelize on parallel computers than direct methods. Different computing platforms,
 to parallelize on parallel computers than direct methods. Different computing platforms,

-sequential and parallel computers, are used for solving sparse linear systems with iterative
-solutions. Nowadays, graphics processing units (GPUs) have become attractive for solving
+sequential and parallel computers, are used to solve sparse linear systems with iterative
+solutions. Nowadays, graphics processing units (GPUs) have become attractive to solve

 these systems, due to their computing power and their ability to compute faster than
 traditional CPUs.
 
 these systems, due to their computing power and their ability to compute faster than
 traditional CPUs.
 


@@ -39,7 +39,7 @@ In Section~\ref{ch12:sec:02}, we describe the general principle of two well-know
 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. Finally, in Section~\ref{ch12:sec:04}, we present the experimental results obtained on a
 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. 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.    
+CPU cluster and on a GPU cluster, to solve large sparse linear systems.    

 
 
 %%--------------------------%%
 
 
 %%--------------------------%%


@@ -81,7 +81,7 @@ linear systems are those called \textit{Krylov subspace methods}~\cite{ch12: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 subspace methods are usually used with preconditioners that allow to improve their
 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}
+convergence. So, in what follows, the CG and GMRES methods are used to solve the left-preconditioned\index{Sparse~linear~system!Preconditioned}

 sparse linear system:
 \begin{equation}
 M^{-1}Ax=M^{-1}b,
 sparse linear system:
 \begin{equation}
 M^{-1}Ax=M^{-1}b,


@@ -94,9 +94,9 @@ where $M$ is the preconditioning matrix.
 %%****************%%
 \subsection{CG method}
 \label{ch12:sec:02.01}
 %%****************%%
 \subsection{CG method}
 \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
+The conjugate gradient method was initially developed by Hestenes and Stiefel in 1952~\cite{ch12:ref2}.
+It is one of the well known iterative method to solve large sparse linear systems. In addition, it
+can be adapted to solve 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
 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


@@ -309,7 +309,7 @@ $V$ to $m$ orthogonal vectors.
 \label{ch12:alg:02}
 \end{algorithm}
 
 \label{ch12:alg:02}
 \end{algorithm}
 

-Algorithm~\ref{ch12:alg:02} shows the main key points of the GMRES method with restarts.
+Algorithm~\ref{ch12:alg:02} shows the 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
 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


@@ -352,7 +352,7 @@ $i$:
 \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
 \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.
+row-wise partitioning (row-by-row decomposition) 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. 
 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. 


@@ -391,10 +391,10 @@ GMRES solvers. The following kernels of CUBLAS (dealing with double floating poi
 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
 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
+develop a kernel for the XPAY operation ($y\leftarrow x+ay$) used  line~$12$ in

 Algorithm~\ref{ch12:alg:01}. In the GMRES solver, we program a kernel for the scalar-vector
 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
+multiplication (lines~$7$ and~$15$ in Algorithm~\ref{ch12:alg:02}), a kernel to
+solve the least-squares problem and a kernel to update the elements  of the solution

 vector $x$.
 
 The least-squares problem in the GMRES method is solved by performing a QR factorization
 vector $x$.
 
 The least-squares problem in the GMRES method is solved by performing a QR factorization


@@ -411,12 +411,12 @@ methods is the sparse matrix-vector multiplication (SpMV)\index{SpMV~multiplicat
 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
 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
+can cause a significant waste of memory space and execution time. In addition, the sparse

 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
 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
+accesses to the global memory, which slows down its performances even more. One of the

 most efficient compressed storage formats\index{Compressed~storage~format} of sparse
 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}.
+matrices on GPUs is the 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
 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


@@ -485,7 +485,7 @@ storage formats. Figure~\ref{ch12:fig:03} shows a reordering of a sparse sub-mat
 
 A GPU cluster\index{GPU~cluster} is a parallel platform with a distributed memory. So, the synchronizations
 and communication data between GPU nodes are carried out by passing messages. However, GPUs can not communicate
 
 A GPU cluster\index{GPU~cluster} is a parallel platform with a distributed memory. So, the synchronizations
 and communication data between GPU nodes are carried out by passing messages. However, GPUs can not communicate

-between them in direct way. Then, CPUs via MPI processes are in charge of the synchronizations within the GPU
+between them in a direct way. Then, CPUs via MPI processes are in charge of the synchronizations within the GPU

 cluster. Consequently, the vector elements to be exchanged must be copied from the GPU memory to the CPU memory
 and vice-versa before and after the synchronization operation between CPUs. We have used the CUBLAS\index{CUBLAS}
 communication subroutines to perform the data transfers between a CPU core and its GPU: \verb+cublasGetVector()+
 cluster. Consequently, the vector elements to be exchanged must be copied from the GPU memory to the CPU memory
 and vice-versa before and after the synchronization operation between CPUs. We have used the CUBLAS\index{CUBLAS}
 communication subroutines to perform the data transfers between a CPU core and its GPU: \verb+cublasGetVector()+


@@ -509,9 +509,9 @@ 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.
 
 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.
 

-Linux cluster version 2.6.39 OS is installed on CPUs. C programming language is used for coding
+Linux cluster version 2.6.39 OS is installed on CPUs. C programming language is used to code

 the parallel algorithms of both methods on the GPU cluster. CUDA version 4.0~\cite{ch12:ref9}
 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
+is used to program GPUs, using CUBLAS library~\cite{ch12:ref6} to deal with vector operations

 in GPUs and, finally, MPI routines of OpenMPI 1.3.3 are used to carry out the communications between
 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.
 in GPUs and, finally, MPI routines of OpenMPI 1.3.3 are used to carry out the communications between
 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.


@@ -573,10 +573,10 @@ over $10$ executions of the same parallel linear solver and for the same input d
 \label{ch12:tab:01}
 \end{table}
 
 \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
+To get more realistic results, we have 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 have chosen six

 symmetric sparse matrices and six nonsymmetric ones from this collection. In Figure~\ref{ch12:fig:05},
 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
+we show the 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.
 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.


@@ -639,11 +639,11 @@ torso3            & $4.242s$           & $2.030s$            & $2.09$        & $
 \end{center}
 \end{table}
 
 \end{center}
 \end{table}
 

-Tables~\ref{ch12:tab:02} and~\ref{ch12:tab:03} shows the performances of the parallel
+Tables~\ref{ch12:tab:02} and~\ref{ch12:tab:03} show the performances of the parallel

 CG and GMRES solvers, respectively, for solving linear systems associated to the sparse
 matrices presented in Tables~\ref{ch12:tab:01}. They allow to compare the performances
 obtained on a cluster of $24$ CPU cores and on a cluster of $12$ GPUs. However, Table~\ref{ch12:tab:02}
 CG and GMRES solvers, respectively, for solving linear systems associated to the sparse
 matrices presented in Tables~\ref{ch12:tab:01}. They allow to compare the performances
 obtained on a cluster of $24$ CPU cores and on a cluster of $12$ GPUs. However, Table~\ref{ch12:tab:02}

-shows only the performances of solving symmetric sparse linear systems, due to the inability
+only shows the performances of solving symmetric sparse linear systems, due to the inability

 of the CG method to solve the nonsymmetric systems. In both tables, the second and third
 columns give, respectively, the execution times in seconds obtained on $24$ CPU cores
 ($Time_{gpu}$) and that obtained on $12$ GPUs ($Time_{gpu}$). Moreover, we take into account
 of the CG method to solve the nonsymmetric systems. In both tables, the second and third
 columns give, respectively, the execution times in seconds obtained on $24$ CPU cores
 ($Time_{gpu}$) and that obtained on $12$ GPUs ($Time_{gpu}$). Moreover, we take into account


@@ -670,18 +670,18 @@ $prec$ is the maximum element, in absolute value, of the residual vector $r^{gpu
 of the solution $x^{gpu}$. Thus, we can see that the solutions obtained on the GPU cluster
 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}$.
 of the solution $x^{gpu}$. Thus, we can see that the solutions obtained on the GPU cluster
 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
+However, we can notice from the relative gains $\tau$ that it is not interesting to use multiple

 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
 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.
+computations over the cluster increase the idle times of GPUs and slow down  the parallel
+computations further.

 
 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 collection~\cite{ch12:ref10}
 
 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 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,
+as an initial matrix to build 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 build its sparse
+sub-matrix. In the first experimental tests, we focused on sparse matrices having a banded structure,
+because they are those arising the most in the majority 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 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
 each MPI process constructs its sub-matrix by performing several copies of an initial sparse matrix chosen
 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


@@ -736,8 +736,8 @@ Tables~\ref{ch12:tab:05} and~\ref{ch12:tab:06} shows the performances of the par
 GMRES solvers, respectively, obtained on a cluster of $24$ CPU cores and on a cluster of $12$
 GPUs. Obviously, we can notice from these tables that solving large sparse linear systems on
 a GPU cluster is more efficient than on a CPU cluster (see relative gains $\tau$). We can also
 GMRES solvers, respectively, obtained on a cluster of $24$ CPU cores and on a cluster of $12$
 GPUs. Obviously, we can notice from these tables that solving large sparse linear systems on
 a GPU cluster is more efficient than on a CPU cluster (see relative gains $\tau$). We can also

-notice that the execution times of the CG method, whether in a CPU cluster or on a GPU cluster,
-are better than those the GMRES method for solving large symmetric linear systems. In fact, the
+notice that the execution times of the CG method, whether in a CPU cluster or in a GPU cluster,
+are better than those of the GMRES method for solving large symmetric linear systems. In fact, the

 CG method is characterized by a better convergence\index{Convergence} rate and a shorter execution
 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.
 CG method is characterized by a better convergence\index{Convergence} rate and a shorter execution
 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.


@@ -810,8 +810,8 @@ on a cluster of 12 GPUs.}
 In this chapter, we have aimed at harnessing the computing power of a
 cluster of GPUs for solving large sparse linear systems. For this, we
 have used two Krylov subspace iterative methods: the CG and GMRES methods.
 In this chapter, we have aimed at harnessing the computing power of a
 cluster of GPUs for solving large sparse linear systems. For this, we
 have used two Krylov subspace iterative methods: the CG and GMRES methods.

-The first method is well-known to its efficiency for solving symmetric
-linear systems and the second one is used, particularly, for solving
+The first method is well-known for its efficiency to solve symmetric
+linear systems and the second one is used, particularly, to solve

 nonsymmetric linear systems. 
 
 We have presented the parallel implementation of both iterative methods
 nonsymmetric linear systems. 
 
 We have presented the parallel implementation of both iterative methods


@@ -825,12 +825,12 @@ algorithms.
 
 In the experimental tests, we have shown that using a GPU cluster is efficient
 for solving linear systems associated to very large sparse matrices. The experimental
 
 In the experimental tests, we have shown that using a GPU cluster is efficient
 for solving linear systems associated to very large sparse matrices. The experimental

-results, obtained in the present chapter, showed that a cluster of $12$ GPUs is
+results, obtained in the present chapter, show that a cluster of $12$ GPUs is

 about $7$ times faster than a cluster of $24$ CPU cores for solving large sparse
 linear systems of $25$ million unknown values. This is due to the GPU ability to
 compute the data-parallel operations faster than the CPUs.
 
 about $7$ times faster than a cluster of $24$ CPU cores for solving large sparse
 linear systems of $25$ million unknown values. This is due to the GPU ability to
 compute the data-parallel operations faster than the CPUs.
 

-As future works, we plan to test the parallel algorithms of CG and GMRES methods, adapted
+In our future works, we plan to test the parallel algorithms of CG and GMRES methods, adapted

 to GPUs, for solving large linear systems associated to sparse matrices of different structures.
 For example, the matrices having large bandwidths, which can lead to many data dependencies
 between the computing nodes and, thus, degrade the performances of both algorithms. So in
 to GPUs, for solving large linear systems associated to sparse matrices of different structures.
 For example, the matrices having large bandwidths, which can lead to many data dependencies
 between the computing nodes and, thus, degrade the performances of both algorithms. So in