X-Git-Url: https://bilbo.iut-bm.univ-fcomte.fr/and/gitweb/kahina_paper1.git/blobdiff_plain/9213a919b9009602014fa0b33a14063f019562dc..fdc11565a5688a74effe8d08b78b4a1853387f3a:/paper.tex?ds=sidebyside diff --git a/paper.tex b/paper.tex index d9b3e5e..694ae43 100644 --- a/paper.tex +++ b/paper.tex @@ -429,74 +429,74 @@ polynomials of 48000. %to $1,000,000$). -\section {A CUDA parallel Ehrlich-Aberth method} -In the following, we describe the parallel implementation of Ehrlich-Aberth method on GPU -for solving high degree polynomials (up to $1,000,000$). First, the hardware and software of the GPUs are presented. Then, the CUDA parallel Ehrlich-Aberth method is presented. - -\subsection{Background on the GPU architecture} -A GPU is viewed as an accelerator for the data-parallel and -intensive arithmetic computations. It draws its computing power -from the parallel nature of its hardware and software -architectures. A GPU is composed of hundreds of Streaming -Processors (SPs) organized in several blocks called Streaming -Multiprocessors (SMs). It also has a memory hierarchy. It has a -private read-write local memory per SP, fast shared memory and -read-only constant and texture caches per SM and a read-write -global memory shared by all its SPs~\cite{NVIDIA10}. - -On a CPU equipped with a GPU, all the data-parallel and intensive -functions of an application running on the CPU are off-loaded onto -the GPU in order to accelerate their computations. A similar -data-parallel function is executed on a GPU as a kernel by -thousands or even millions of parallel threads, grouped together -as a grid of thread blocks. Therefore, each SM of the GPU executes -one or more thread blocks in SIMD fashion (Single Instruction, -Multiple Data) and in turn each SP of a GPU SM runs one or more -threads within a block in SIMT fashion (Single Instruction, -Multiple threads). Indeed at any given clock cycle, the threads -execute the same instruction of a kernel, but each of them -operates on different data. - GPUs only work on data filled in their -global memories and the final results of their kernel executions -must be communicated to their CPUs. Hence, the data must be -transferred in and out of the GPU. However, the speed of memory -copy between the GPU and the CPU is slower than the memory -bandwidths of the GPU memories and, thus, it dramatically affects -the performances of GPU computations. Accordingly, it is necessary -to limit as much as possible, data transfers between the GPU and its CPU during the -computations. -\subsection{Background on the CUDA Programming Model} - -The CUDA programming model is similar in style to a single program -multiple-data (SPMD) software model. The GPU is viewed as a -coprocessor that executes data-parallel kernel functions. CUDA -provides three key abstractions, a hierarchy of thread groups, -shared memories, and barrier synchronization. Threads have a three -level hierarchy. A grid is a set of thread blocks that execute a -kernel function. Each grid consists of blocks of threads. Each -block is composed of hundreds of threads. Threads within one block -can share data using shared memory and can be synchronized at a -barrier. All threads within a block are executed concurrently on a -multithreaded architecture.The programmer specifies the number of -threads per block, and the number of blocks per grid. A thread in -the CUDA programming language is much lighter weight than a thread -in traditional operating systems. A thread in CUDA typically -processes one data element at a time. The CUDA programming model -has two shared read-write memory spaces, the shared memory space -and the global memory space. The shared memory is local to a block -and the global memory space is accessible by all blocks. CUDA also -provides two read-only memory spaces, the constant space and the -texture space, which reside in external DRAM, and are accessed via -read-only caches. - -\section{ The implementation of Ehrlich-Aberth method on GPU} +%% \section {A CUDA parallel Ehrlich-Aberth method} +%% In the following, we describe the parallel implementation of Ehrlich-Aberth method on GPU +%% for solving high degree polynomials (up to $1,000,000$). First, the hardware and software of the GPUs are presented. Then, the CUDA parallel Ehrlich-Aberth method is presented. + +%% \subsection{Background on the GPU architecture} +%% A GPU is viewed as an accelerator for the data-parallel and +%% intensive arithmetic computations. It draws its computing power +%% from the parallel nature of its hardware and software +%% architectures. A GPU is composed of hundreds of Streaming +%% Processors (SPs) organized in several blocks called Streaming +%% Multiprocessors (SMs). It also has a memory hierarchy. It has a +%% private read-write local memory per SP, fast shared memory and +%% read-only constant and texture caches per SM and a read-write +%% global memory shared by all its SPs~\cite{NVIDIA10}. + +%% On a CPU equipped with a GPU, all the data-parallel and intensive +%% functions of an application running on the CPU are off-loaded onto +%% the GPU in order to accelerate their computations. A similar +%% data-parallel function is executed on a GPU as a kernel by +%% thousands or even millions of parallel threads, grouped together +%% as a grid of thread blocks. Therefore, each SM of the GPU executes +%% one or more thread blocks in SIMD fashion (Single Instruction, +%% Multiple Data) and in turn each SP of a GPU SM runs one or more +%% threads within a block in SIMT fashion (Single Instruction, +%% Multiple threads). Indeed at any given clock cycle, the threads +%% execute the same instruction of a kernel, but each of them +%% operates on different data. +%% GPUs only work on data filled in their +%% global memories and the final results of their kernel executions +%% must be communicated to their CPUs. Hence, the data must be +%% transferred in and out of the GPU. However, the speed of memory +%% copy between the GPU and the CPU is slower than the memory +%% bandwidths of the GPU memories and, thus, it dramatically affects +%% the performances of GPU computations. Accordingly, it is necessary +%% to limit as much as possible, data transfers between the GPU and its CPU during the +%% computations. +%% \subsection{Background on the CUDA Programming Model} + +%% The CUDA programming model is similar in style to a single program +%% multiple-data (SPMD) software model. The GPU is viewed as a +%% coprocessor that executes data-parallel kernel functions. CUDA +%% provides three key abstractions, a hierarchy of thread groups, +%% shared memories, and barrier synchronization. Threads have a three +%% level hierarchy. A grid is a set of thread blocks that execute a +%% kernel function. Each grid consists of blocks of threads. Each +%% block is composed of hundreds of threads. Threads within one block +%% can share data using shared memory and can be synchronized at a +%% barrier. All threads within a block are executed concurrently on a +%% multithreaded architecture.The programmer specifies the number of +%% threads per block, and the number of blocks per grid. A thread in +%% the CUDA programming language is much lighter weight than a thread +%% in traditional operating systems. A thread in CUDA typically +%% processes one data element at a time. The CUDA programming model +%% has two shared read-write memory spaces, the shared memory space +%% and the global memory space. The shared memory is local to a block +%% and the global memory space is accessible by all blocks. CUDA also +%% provides two read-only memory spaces, the constant space and the +%% texture space, which reside in external DRAM, and are accessed via +%% read-only caches. + +\section{ Implementation of Ehrlich-Aberth method on GPU} \label{sec5} %%\subsection{A CUDA implementation of the Aberth's method } %%\subsection{A GPU implementation of the Aberth's method } -\subsection{A sequential Ehrlich-Aberth algorithm} +\subsection{Sequential Ehrlich-Aberth algorithm} The main steps of Ehrlich-Aberth method are shown in Algorithm.~\ref{alg1-seq} : %\LinesNumbered \begin{algorithm}[H] @@ -504,8 +504,10 @@ The main steps of Ehrlich-Aberth method are shown in Algorithm.~\ref{alg1-seq} : \caption{A sequential algorithm to find roots with the Ehrlich-Aberth method} -\KwIn{$Z^{0}$(Initial root's vector),$\varepsilon$ (error tolerance threshold), P(Polynomial to solve),$\Delta z_{max}$ (maximum value of stop condition),k (number of iteration),n(Polynomial's degrees)} -\KwOut {Z (The solution root's vector),ZPrec (the previous solution root's vector)} +\KwIn{$Z^{0}$ (Initial root's vector), $\varepsilon$ (error tolerance + threshold), P (Polynomial to solve), $\Delta z_{max}$ (maximum value + of stop condition), k (number of iteration), n (Polynomial's degrees)} +\KwOut {Z (The solution root's vector), ZPrec (the previous solution root's vector)} \BlankLine @@ -546,17 +548,24 @@ With the Gauss-Seidel iteration, we have: EAGS: z^{k+1}_{i}=\frac{p(z^{k}_{i})}{p'(z^{k}_{i})-p(z^{k}_{i})(\sum^{i-1}_{j=1}\frac{1}{z^{k}_{i}-z^{k+1}_{j}}+\sum^{n}_{j=i+1}\frac{1}{z^{k}_{i}-z^{k}_{j}})}, i=1,...,n. \end{equation} %%Here a finiched my revision %% -Using Equation.~\ref{eq:Aberth-H-GS} to update the vector solution \textit{Z}, we expect the Gauss-Seidel iteration to converge more quickly because, just as its ancestor (for solving linear systems of equations), it uses the most fresh computed roots $z^{k+1}_{i}$. +Using Equation.~\ref{eq:Aberth-H-GS} to update the vector solution +\textit{Z}, we expect the Gauss-Seidel iteration to converge more +quickly because, just as any Jacobi algorithm (for solving linear systems of equations), it uses the most fresh computed roots $z^{k+1}_{i}$. The $4^{th}$ step of the algorithm checks the convergence condition using Equation.~\ref{eq:Aberth-Conv-Cond}. Both steps 3 and 4 use 1 thread to compute all the $n$ roots on CPU, which is very harmful for performance in case of the large degree polynomials. -\subsection{A Parallel implementation with CUDA } +\subsection{Parallel implementation with CUDA } On the CPU, both steps 3 and 4 contain the loop \verb=for= and a single thread executes all the instructions in the loop $n$ times. In this subsection, we explain how the GPU architecture can compute this loop and reduce the execution time. -In the GPU, the schduler assigns the execution of this loop to a group of threads organised as a grid of blocks with block containing a number of threads. All threads within a block are executed concurrently in parallel. The instructions run on the GPU are grouped in special function called kernels. It's up to the programmer, to describe the execution context, that is the size of the Grid, the number of blocks and the number of threads per block upon the call of a given kernel, according to a special syntax defined by CUDA. +In the GPU, the scheduler assigns the execution of this loop to a +group of threads organised as a grid of blocks with block containing a +number of threads. All threads within a block are executed +concurrently in parallel. The instructions run on the GPU are grouped +in special function called kernels. With CUDA, a programmer must +describe the kernel execution context: the size of the Grid, the number of blocks and the number of threads per block. -In CUDA programming, all the instructions of the \verb=for= loop are executed by the GPU as a kernel. A kernel is a function written in CUDA and defined by the \verb=__global__= qualifier added before a usual \verb=C= function, which instructs the compiler to generate appropriate code to pass it to the CUDA runtime in order to be executed on the GPU. +%In CUDA programming, all the instructions of the \verb=for= loop are executed by the GPU as a kernel. A kernel is a function written in CUDA and defined by the \verb=__global__= qualifier added before a usual \verb=C= function, which instructs the compiler to generate appropriate code to pass it to the CUDA runtime in order to be executed on the GPU. Algorithm~\ref{alg2-cuda} shows a sketch of the Ehrlich-Aberth algorithm using CUDA. @@ -565,13 +574,13 @@ Algorithm~\ref{alg2-cuda} shows a sketch of the Ehrlich-Aberth algorithm using C %\LinesNumbered \caption{CUDA Algorithm to find roots with the Ehrlich-Aberth method} -\KwIn{$Z^{0}$(Initial root's vector),$\varepsilon$ (error tolerance threshold), P(Polynomial to solve), $\Delta z_{max}$ (maximum value of stop condition)} +\KwIn{$Z^{0}$ (Initial root's vector), $\varepsilon$ (error tolerance threshold), P(Polynomial to solve), $\Delta z_{max}$ (maximum value of stop condition)} -\KwOut {Z(The solution root's vector)} +\KwOut {Z (The solution root's vector)} \BlankLine -Initialization of the coeffcients of the polynomial to solve\; +Initialization of the coefficients of the polynomial to solve\; Initialization of the solution vector $Z^{0}$\; Allocate and copy initial data to the GPU global memory\; k=0\; @@ -588,7 +597,11 @@ $kernel\_testConverge(\Delta z_{max},d\_Z^{k},d\_Z^{k-1})$\; After the initialisation step, all data of the root finding problem to be solved must be copied from the CPU memory to the GPU global memory, because the GPUs only access data already present in their memories. Next, all the data-parallel arithmetic operations inside the main loop \verb=(do ... while(...))= are executed as kernels by the GPU. The first kernel named \textit{save} in line 6 of Algorithm~\ref{alg2-cuda} consists in saving the vector of polynomial's root found at the previous time-step in GPU memory, in order to check the convergence of the roots after each iteration (line 8, Algorithm~\ref{alg2-cuda}). -The second kernel executes the iterative function $H$ and updates $z^{k}$, according to Algorithm~\ref{alg3-update}. We notice that the update kernel is called in two forms, separated with the value of \emph{R} which determines the radius beyond which we apply the logarithm computation of the power of a complex. +The second kernel executes the iterative function $H$ and updates +$z^{k}$, according to Algorithm~\ref{alg3-update}. We notice that the +update kernel is called in two forms, separated with the value of +\emph{R} which determines the radius beyond which we apply the +exponential logarithm algorithm. \begin{algorithm}[H] \label{alg3-update} @@ -598,18 +611,38 @@ The second kernel executes the iterative function $H$ and updates $z^{k}$, accor \eIf{$(\left|Z^{(k)}\right|<= R)$}{ $kernel\_update(d\_z^{k})$\;} { -$kernel\_update\_Log(d\_z^{k})$\; +$kernel\_update\_ExpoLog(d\_z^{k})$\; } \end{algorithm} -The first form executes formula \ref{eq:SimplePolynome} if the modulus of the current complex is less than the a certain value called the radius i.e. ($ |z^{k}_{i}|<= R$), else the kernel executes formulas (Eq.~\ref{deflncomplex},Eq.~\ref{defexpcomplex}). The radius $R$ is evaluated as : +The first form executes formula \ref{eq:SimplePolynome} if the modulus +of the current complex is less than the a certain value called the +radius i.e. ($ |z^{k}_{i}|<= R$), else the kernel executes the EA.EL +function Eq.~\ref{Log_H2} +(with Eq.~\ref{deflncomplex}, Eq.~\ref{defexpcomplex}). The radius $R$ is evaluated as : $$R = \exp( \log(DBL\_MAX) / (2*n) )$$ where $DBL\_MAX$ stands for the maximum representable double value. -The last kernel verifies the convergence of the roots after each update of $Z^{(k)}$, according to formula. We used the functions of the CUBLAS Library (CUDA Basic Linear Algebra Subroutines) to implement this kernel. +The last kernel checks the convergence of the roots after each update +of $Z^{(k)}$, according to formula Eq.~\ref{eq:Aberth-Conv-Cond}. We used the functions of the CUBLAS Library (CUDA Basic Linear Algebra Subroutines) to implement this kernel. + +The kernel terminates its computations when all the roots have +converged. Many important remarks should be noticed. First, as blocks +of threads are scheduled automatically by the GPU, we have absolutely +no control on the order of the blocks. Consequently, our algorithm is +executed more or less in an asynchronous iterations way, where blocks +of roots are updated in a non deterministic way. As the Durand-Kerner +method has been proved to convergence with asynchronous iterations, we +think it is similar with the Ehrlich-Aberth method, but we did not try +to prove this in that paper. Another consequence of that, is that +several executions of our algorithm with the same polynomials do no +give necessarily the same result with the same number of iterations +(even if the variation is not very significant). + + + + -The kernels terminate it computations when all the roots converge. Finally, the solution of the root finding problem is copied back from GPU global memory to CPU memory. We use the communication functions of CUDA for the memory allocation in the GPU \verb=(cudaMalloc())= and for data transfers from the CPU memory to the GPU memory \verb=(cudaMemcpyHostToDevice)= -or from GPU memory to CPU memory \verb=(cudaMemcpyDeviceToHost))=. %%HIER END MY REVISIONS (SIDER) \section{Experimental study} \label{sec6}